SciPyのF分布ppfが抱える数値精度問題

概要 scipy.stats.f.ppf(F分布のパーセント点関数)が、特定の条件下で不正確な結果を返す可能性があるという問題についての報告です。具体的には、確率が1に非常に近い値(裾の領域)と大きな自由度の組み合わせにおいて、ppf関数の結果と、累積分布関数(cdf)を最適化して得られる結果との間に大きな乖離が見られます。 この問題は2021年頃に発見されましたが、最近のバージョンでも同様の現象が確認されたため、改めて報告するものです。 なお、本件は、既に、https://github.com/scipy/scipy/issues/20835 で報告されており、恐らく、scipy 1,17で対策はマージされると思います。 ENH: special: boostify F distribution functions 問題の詳細 scipy.stats.f.ppf に欠陥がある可能性が疑われます。最適化関数と累積分布関数(cdf)の組み合わせで得られる結果が、ppf関数の戻り値と一致しません。 再現コード 以下のコードは、sps.f.ppf を直接呼び出した場合と、spo.brentq を使って sps.f.cdf から逆関数を求めた場合の結果を比較します。 import scipy.stats as sps import scipy.optimize as spo # ppf関数を直接呼び出す print(sps.f.ppf(0.99999995, 1, 50000, loc=0, scale=1)) # cdfと最適化関数brentqを使って逆関数を求める print(spo.brentq(lambda x: sps.f.cdf(x, 1, 50000) - 0.99999995, 1, 1e10)) コードの出力 3333.8385803475894 29.725915444860586 出力結果からわかるように、2つの計算結果は大きく異なっています。 エラーメッセージ この問題では、エラーメッセージは出力されません。コードは正常に終了しますが、得られる結果の信頼性に疑問があります。 原因の考察 このPythonコードと出力は、F分布のパーセント点関数(ppf)と、それを最適化関数(brentq)を使って求めるアプローチとの間の数値的な不安定性を示唆しています。 sps.f.ppf(0.99999995, 1, 50000) これはF分布のパーセント点関数(ppf)を直接呼び出し、累積確率が0.99999995となるF値を求めます。この計算では、F値が 3333.83… という非常に大きな値になっています。 spo.brentq(lambda x: sps.f.cdf(x, 1, 50000) - 0.99999995, 1, 1e10) こちらはscipy.optimizeモジュールのbrentq関数を使い、同じ問題を数値的に解いています。brentqは、与えられた関数の値が0になる点(根)を見つけるための堅牢なアルゴリズムです。ここでは、「sps.f.cdf(x, 1, 50000) - 0.99999995」という式が0になるx、つまり累積確率が0.99999995になるF値を探索しています。 この計算では、F値が 29.72… という、ppfの直接呼び出しとは全く異なる、より妥当と思われる値が得られます。 この乖離は、scipy.stats.f.ppfが特定の条件下で正確な答えを返さない可能性があることを示しています。 F分布のような統計分布の逆関数(ppf)を、極端な確率値(この例では1に非常に近い0.99999995)で計算する場合、数値的な精度の問題が発生しやすくなります。特に自由度が大きい場合、分布の裾(テール)の部分が非常に平坦になるため、計算過程でのわずかな丸め誤差などが、結果として得られるF値の大きな違いにつながることがあります。 一方、brentqは指定された範囲内で根を探索する、より頑健な数値計算アルゴリズムです。関数の振る舞いをより正確に追跡できるため、このようなケースではより信頼性の高い結果を生成すると考えられます。 問題の可視化 この問題は、以下のコードでscipy.special.fdtri(F分布の逆関数)の挙動をプロットすることで、より明確に可視化できます。自由度d2が50000と大きい場合に、確率pが0に近づく(つまり累積確率が1に近づく)領域で、関数の値が急激に変動していることがわかります。 import numpy as np from scipy import special import matplotlib.pyplot as plt d1, d2, p = 1, 50000, np.logspace(-16, -1) # fdtriは1-pを引数に取るため、pが0に近づくと累積確率が1に近づく fdtri = special.fdtri(d1, d2, 1-p) plt.semilogx(p, fdtri) plt.xlabel("p (1 - Cumulative Probability)") plt.ylabel("F-value (fdtri)") plt.title("fdtri behavior for extreme probabilities (d1=1, d2=50000)") plt.grid(True) plt.show() ...

8月 19, 2025 · 2 分 · 279 文字 · Me

Record of Signate War

“採血データを使った心不全予測"のコンペに参加したので、解法を示します。まず、基本的に、2つの分類の予測問題は元々慣れていたので基本的なフレームとしては以下の通りです。 探索的データ分析 ベースラインモデルの作成 予測モデルの作成 探索的データ分析 まず、最初にデータを一応概観してみました。 発売中の技術同人誌 “Pythonによる探索的データ分析クックブック“でも触れている、ydata-profilingを使用しています。 import os import sys import pandas as pd import polars as pl import pyarrow as pa import numpy as np import matplotlib.pyplot as plt import seaborn as sns from ydata_profiling import ProfileReport train_df = pd.read_csv("../data/train.csv") test_df = pd.read_csv("../data/test.csv") profile = ProfileReport(train_df, title="Heart Failure Report") profile.to_file("../profile/heart_failure_report.html") ベースラインモデルの作成 基本線となるベースラインモデルを作成しました。ベースラインモデルは文字通りベースラインモデルなので、複雑なものは避けるのがセオリーです。今回は二つの分類をするタイプなので、一般化線形回帰でロジスティック回帰に持ち込むのを基本としました。また、stepwiseなどの容易さから、一旦、Rでモデリングを進めました。 コードとしては以下のシンプルきわまるものです。 require(dplyr) require(readr) require(ggplot2) require(pROC) train.df <- read.csv("../data/train.csv") test.df <- read.csv("../data/test.csv") train.df <- train.df %>% mutate( anaemia = as.factor(anaemia), diabetes = as.factor(diabetes), high_blood_pressure = as.factor(high_blood_pressure), sex = as.factor(sex), smoking = as.factor(smoking), target = as.factor(target) ) require(MASS) model <- glm(formula = target ~ age + anaemia + creatinine_phosphokinase + diabetes + ejection_fraction + high_blood_pressure + platelets + serum_creatinine + serum_sodium + sex + smoking + time, data = train.df, family = binomial) model.opt <- stepAIC(model) pred.df <- predict(model.opt, train.df, type = "response") roc_curve <- roc(train.df$target, pred.df) auc_value <- auc(roc_curve) ar_value <- (auc_value - 0.5) * 2 print(paste("AR値:", ar_value)) cat("Length of Cumulative_Percentage:", length(seq(0, 1, length.out = nrow(train.df))), "\n") cat("Length of Predicted_Positive:", length(sort(pred.df, decreasing = TRUE)), "\n") cat("Length of Actual_Positive:", length(sort(as.numeric(train.df$target), decreasing = TRUE)), "\n") sorted_predictions <- sort(pred.df, decreasing = TRUE) # 予測確率を降順に sorted_targets <- sort(as.numeric(as.character(train.df$target)), decreasing = TRUE) # 実ターゲットを降順に n <- min(length(sorted_predictions), length(sorted_targets)) cumulative_data <- data.frame( Cumulative_Percentage = seq(0, 1, length.out = n), Predicted_Positive = cumsum(sorted_predictions[1:n]), Actual_Positive = cumsum(sorted_targets[1:n]) ) ggplot(cumulative_data, aes(x = Cumulative_Percentage)) + geom_line(aes(y = Predicted_Positive, color = "モデル予測")) + geom_line(aes(y = Actual_Positive, color = "実データ")) + scale_y_continuous(name = "累積正例率", limits = c(0, max(cumulative_data$Predicted_Positive, cumulative_data$Actual_Positive))) + scale_x_continuous(name = "累積パーセンテージ", limits = c(0, 1)) + labs(title = "CAP図") + theme_minimal() 最終的には以下のCAP図になりました、 ...

5月 2, 2025 · 5 分 · 1037 文字 · Me

RinnaStableDiffusion

Rinnaから日本語対応のStable Diffusionが出たのでをGoogle Colab上で使ってみました。 コードとしては以下のような形です。bashのコードはJupyterから投げます。 pip install gradio try: from japanese_stable_diffusion import JapaneseStableDiffusionPipeline except: res = subprocess.run(['pip', 'install', 'git+https://github.com/rinnakk/japanese-stable-diffusion'], stdout=subprocess.PIPE).stdout.decode('utf-8') print(res) from japanese_stable_diffusion import JapaneseStableDiffusionPipeline import torch from torch import autocast from diffusers import LMSDiscreteScheduler from PIL import Image from IPython import display import gradio as gr def make_grid_from_pils(pil_images): w, h = pil_images[0].size grid_img = Image.new("RGB", ((len(pil_images)) * w, h)) for idx, image in enumerate(pil_images): grid_img.paste(image, (idx * w, 0)) return grid_img from huggingface_hub import notebook_login notebook_login() model_id = "rinna/japanese-stable-diffusion" device = "cuda" if torch.cuda.is_available() else "cpu" # Use the K-LMS scheduler here instead scheduler = LMSDiscreteScheduler( beta_start=0.00085, beta_end=0.012, beta_schedule="scaled_linear", num_train_timesteps=1000 ) pipe = JapaneseStableDiffusionPipeline.from_pretrained( pretrained_model_name_or_path=model_id, scheduler=scheduler, torch_dtype=torch.float16, use_auth_token=True ).to(device) #@markdown ###**Inference Setting:** # the number of output images. If you encounter Out Of Memory error, decrease this number. n_samples = 1 #@param{type: 'integer'} # `classifier-free guidance scale` adjusts how much the image will be like your prompt. Higher values keep your image closer to your prompt. guidance_scale = 7.5 #@param {type:"number"} # How many steps to spend generating (diffusing) your image. steps = 50 #@param{type: 'integer'} # The width of the generated image. width = 512 #@param{type: 'integer'} # The height of the generated image. height = 512 #@param{type: 'integer'} # The seed used to generate your image. Enable to manually set a seed. seed = 'random' #@param{type: 'string'} import torch from torch import autocast from diffusers import LMSDiscreteScheduler from japanese_stable_diffusion import JapaneseStableDiffusionPipeline model_id = "rinna/japanese-stable-diffusion" device = "cuda" # Use the K-LMS scheduler here instead scheduler = LMSDiscreteScheduler(beta_start=0.00085, beta_end=0.012, beta_schedule="scaled_linear", num_train_timesteps=1000) pipe = JapaneseStableDiffusionPipeline.from_pretrained(model_id, scheduler=scheduler, use_auth_token=True) pipe = pipe.to(device) prompt = "富士山をバックに二大スーパーロボットががっちりと握手" with autocast("cuda"): image = pipe(prompt, guidance_scale=7.5)["sample"][0] image.save("output.png") image “富士山をバックに二大スーパーロボットががっちりと握手"から画像を作成し、以下のような画像になります。 ...

9月 9, 2022 · 2 分 · 272 文字 · Me

StableDiffusion

Stable DiffusionをGoogle Colab上で使ってみました。 コードとしては以下のような形です。bashのコードはJupyterから投げます。 pip install diffusers==0.2.4 pip install transformers scipy ftfy pip install "ipywidgets>=7,<8" from google.colab import output output.enable_custom_widget_manager() from huggingface_hub import notebook_login notebook_login() import torch from diffusers import StableDiffusionPipeline # make sure you're logged in with `huggingface-cli login` pipe = StableDiffusionPipeline.from_pretrained("CompVis/stable-diffusion-v1-4", revision="fp16", torch_dtype=torch.float16, use_auth_token=True) pipe = pipe.to("cuda") from torch import autocast prompt = "a galaxy far from earth" with autocast("cuda"): image = pipe(prompt)["sample"][0] # image here is in [PIL format](https://pillow.readthedocs.io/en/stable/) # Now to display an image you can do either save it such as: image.save(f"galaxy_far_from_earth.png") # or if you're in a google colab you can directly display it with image “a galaxy far from earth"から画像を作成し、以下のような画像になります。 ...

9月 1, 2022 · 1 分 · 115 文字 · Me

Transformersによる文書の分類

Hugging Face Transformersを使ってネガポジを判定するモデルを作ってみました。 query title label negaposi この映画は本当に面白い 0 みたいな形で教師を作り、それを投入して学習させました。 東北大学の日本語 BERT モデルを事前学習モデルとし、 それをSequence Classificationさせました。 モデリング自体は、Google Colaboratoryを用いて実行しました。 学習 !pip install transformers[ja]==4.3.3 torch==1.9 sentencepiece==0.1.91 from google.colab import drive import pandas as pd from sklearn.model_selection import train_test_split from transformers import BertJapaneseTokenizer, BertForSequenceClassification, BertForMaskedLM, pipeline, Trainer, TrainingArguments import torch drive.mount('/content/drive') training_data = pd.read_csv('/content/drive/MyDrive/Texts/negaposi-sentence.csv') training_data.head() print(len(training_data["query"].unique())) training_data[["title", "label"]].groupby("label").count() train_queries, val_queries, train_docs, val_docs, train_labels, val_labels = train_test_split( training_data["query"].tolist(), training_data["title"].tolist(), training_data["label"].tolist(), test_size=.5 ) model_name = 'cl-tohoku/bert-base-japanese-whole-word-masking' tokenizer = BertJapaneseTokenizer.from_pretrained(model_name) train_encodings = tokenizer(train_queries, train_docs, truncation=True, padding='max_length', max_length=128) val_encodings = tokenizer(val_queries, val_docs, truncation=True, padding='max_length', max_length=128) model_name = 'cl-tohoku/bert-base-japanese-whole-word-masking' tokenizer = BertJapaneseTokenizer.from_pretrained(model_name) train_encodings = tokenizer(train_queries, train_docs, truncation=True, padding='max_length', max_length=128) val_encodings = tokenizer(val_queries, val_docs, truncation=True, padding='max_length', max_length=128) model = BertForSequenceClassification.from_pretrained(model_name, num_labels = 2) for param in model.base_model.parameters(): param.requires_grad = False training_args = TrainingArguments( logging_steps=10, output_dir='models', evaluation_strategy="epoch", num_train_epochs=2000, per_device_train_batch_size=16, per_device_eval_batch_size=64, warmup_steps=500, weight_decay=0.01, save_total_limit=1, ) trainer = Trainer( model=model, args=training_args, train_dataset=train_dataset, eval_dataset=val_dataset ) trainer.train() trainer.save_model(output_dir='/content/drive/MyDrive/Models/sentiment-mining4') 推論 !pip install transformers[ja]==4.3.3 torch==1.9 sentencepiece==0.1.91 from google.colab import drive import pandas as pd from sklearn.model_selection import train_test_split from transformers import BertJapaneseTokenizer, BertForSequenceClassification, BertForMaskedLM, pipeline, Trainer, TrainingArguments import torch drive.mount('/content/drive') model_name = 'cl-tohoku/bert-base-japanese-whole-word-masking' tokenizer = BertJapaneseTokenizer.from_pretrained(model_name) model = BertForSequenceClassification.from_pretrained('/content/drive/MyDrive/Models/sentiment-mining4') nlp = pipeline("sentiment-analysis",model=model,tokenizer=tokenizer) nlp("この本は興味深い") nlp(“この本は興味深い”) ...

12月 11, 2021 · 1 分 · 212 文字 · Me