確率密度関数に基づく乱数を生成

Question

前提・実現したいこと

自分で定義した確率密度関数に基づく乱数を生成

発生している問題・エラーメッセージ

そもそもアプローチの仕方がわからない．．
scipyのrvsもここでは使えなさそうだった．

該当のソースコード

import numpy as np
import random
lamda=1.5
xi=0.4
x=random.uniform(1.01,1.5)#違うというかいらない気がする
pdf=1/np.sqrt(2*np.pi)/xi/(x-1)*np.exp(-((np.log(x-1)-lamda)/xi)**2/2)
#このあとがどうしたらいいかわからない．
#random(pdf)???

試したこと

scipy.stats.rvs_ratio_uniforms：ドキュメンテーション読んで使えそうだったのにimport時に下記のようなエラー
import scipy.stats.rvs_ratio_uniforms as rvs

ModuleNotFoundError: No module named 'scipy.stats.rvs_ratio_uniforms'

補足情報（FW/ツールのバージョンなど）

環境はanacondaです．
pdf(x)とxの範囲は定義できています．xの定義域は，1<=x<=無限．x=1.2程度でピークの非対称の分布．
分布は対数正規分布をベースにしたものですが，平均・標準偏差を改良してあるため，一般的なscipyのモジュールなどは使えません．
正規分布の例などあれば，自分で直していけると思います．
あるいは，前述のモジュールを使えるようになる方法を教えてください．
お願いします．

Comments

[hayataka2049]

import時に出たエラーを記載してください。

[hayataka2049]

どんな分布をやる想定なのかも書いていただけると回答しやすいです。

[aochan01]

import scipy.stats.rvs_ratio_uniforms as rvs

ModuleNotFoundError: No module named 'scipy.stats.rvs_ratio_uniforms'
です

[hayataka2049]

他の人も見やすいように質問を編集してそちらにも書いておいてください。

[aochan01]

分布は対数正規分布を改良したものです

Answer 1

解決済みですが、面白そうだったのでスクラッチで実装してみました。

import numpy as np
lamda=1.5
xi=0.4

x = np.arange(1.01,15,0.001)
pdf=1/np.sqrt(2*np.pi)/xi/(x-1)*np.exp(-((np.log(x-1)-lamda)/xi)**2/2)

y_lower = np.min(pdf)
y_upper = np.max(pdf)

trial = 10000
res = []
for i in range(trial) :
    x = np.random.uniform(1.01,15)
    y = np.random.uniform(y_lower,y_upper)
    pdf = 1/np.sqrt(2*np.pi)/xi/(x-1)*np.exp(-((np.log(x-1)-lamda)/xi)**2/2)
    
    if y<=pdf:
        res.append(x)

モンテカルロ法で円周率を計算するコードから着想を得て、確率分布の理論値よりも上か下かを判定するという方法で乱数を生成しました。
結果であるresのヒストグラムを作成したところ、自己解決にあるものに近似したものが出力しました。

Answer 2

importは結局できませんでしたが，random.lognormvariateの数値をいじっくていたら，欲しい値が得られました．協力ありがとうございました．

Answer 3

回答ではありませんが参考までに・・・
pythonを使って、確率密度関数(正規分布、標準正規分布)を覚えちゃう！

＃“確率密度関数 python”で、ググればすぐ出ますよｄ＾＾

Comments

[aochan01]

正規分布やポアソン分布など一般的な分布であれば，scipy.stats.normとかで余裕なんですが，自分で定義した分布なので，この方法が使えなくて困っている次第です．．

Answer 4

以下のQiitaの記事が回答になっていると思います。
Scipyもrvsも使えます。

Qiita Pythonで任意の確率密度関数からサンプリングを行う方法

Answer 5

とりあえず、rvs_ratio_uniformsを使いたい場合はscipy 1.2.0以降が要求されます。環境が古くないかを確認してください。

SciPy 1.2.0 Release Notes — SciPy v1.2.0 Reference Guide

また、以下の書き方でimportしてください。

from scipy.stats import rvs_ratio_uniforms as rvs

Comments

[aochan01]

更新してimportを変えてみましたが，依然エラーがでます．
ImportError: cannot import name 'rvs_ratio_uniforms' from 'scipy.stats' (C:\ProgramData\Anaconda3\lib\site-packages\scipy\stats\__init__.py)

[hayataka2049]

import scipy
して
scipy.version.version
を確認してください。1.2.0になりますか？

[aochan01]

自分，モジュールの入れ方や確認法がよくわからないんですけど，anaconda navigaterによれば，scipyは1.2.0と出ていました．
確率密度関数の定義式から，乱数を生成する方法（参考のページ等）を教えていただければ，なんとかなりそうなんですが．