Category: Math

[C言語]素因数分解

素因数分解とは自然数を素数の掛け算で表すこと
素因数分解の対象は自然数であること(0を除く)
素数とは、1とその数以外に約数を持たないこと

### 素因数分解のやり方？
対象の自然数をNとする
1. Nを2で割る
-> 2で割り切れた場合: 割り切れた値を更に2で割る
-> 2で割り切れない場合: 下に行く
↓
2.3で割る(※2+1)
-> 3で割り切れた場合: 割り切れた値を更に3で割る
-> 3で割り切れない場合: 下に行く
↓
3.5で割る(※2+2)
-> 5で割り切れた場合: 割り切れた値を更に5で割る
-> 5で割り切れない場合: 7 下に行く

2から始めて、これをN回繰り返す　であってる？？
割る数は2から3へは+1だが、3以降は+2の方が無駄な計算がなくなる
プログラムで書きたい

#include 

int main(void){
	int prime = 27; // 素因数
	int assemble[prime]; // 割り切れる数
	int i, j, k;

	j = 0;
	for(i=2; i < prime; i++){
		while(prime % i == 0){
			prime = prime / i;
			assemble[j] = i;
			j++;
		}
	}

	for(k=0; k < j; k++){
		printf("%d ", assemble[k]);
	}
	printf("\n");

	return 0;
}

素因数が27の時
$ ./main
3 3 3
上手くいってる。

素因数が21の時
$ ./main
3

あれ、何でだ。。
あ、forループでprimeが割られた数を代入するからおかしくなるんだ。
forループの変数を変えて再度計算します。

#include 

int main(void){
	int prime = 60; // 素因数
	int cal = 60;
	int assemble[prime]; // 割り切れる数
	int i, j, k;

	j = 0;
	for(i=2; i < cal; i++){
		while(prime % i == 0){
			prime = prime / i;
			assemble[j] = i;
			j++;
		}
	}

	for(k=0; k < j; k++){
		printf("%d ", assemble[k]);
	}
	printf("\n");

	return 0;
}

素因数が60の時
$ ./main
2 2 3 5

素因数が115の時
$ ./main
5 23

わお、中々素晴らしい。
forループの時に、最初だけi++として、2回目以降をi+２とするにはどう書けばいいんだろう。
1回目の割り算を切り離して書いて、2回目以降をforループで書くのだろうか。

クロスエントロピー誤差

目的関数J(Θ)を最小にするには、
J(Θ)= 1/2m*mΣi=1(pxi – yi)^2

ただし、0≦pxi≦1でyiは0もしくは1の為、
L(Θ)= – mΣi=1(yilog(pxi)+(1-yi)log(1-pxi))を求めた方がより簡便な計算になる。
これをクロスエントロピー誤差関数という。

def cross_entropy_error(y,t):
	delta = le-7
	return -np.sum(t*np.log(y+delta))

うん、ここは後々復習が必要そうだな。

最小二乗法

まず、散布図を作成します。

import matplotlib.pyplot as plt 

x = [1,2,3,4,5]
y = [1571,1630,1681,1603,1623]

plt.scatter(x, y)

合計距離の最小値を求める
D = 5Σ[t=1]*{yl – (w0 + w1xl)}^2

実装する

import numpy as np
import matplotlib.pyplot as plt 

x = np.array([1,2,3,4,5])
y = np.array([1571,1630,1681,1603,1623])

def reg1dim(x, y):
	a = np.dot(x, y)/ (x ** 2).sum()
	return a 

a = reg1dim(x, y)

plt.scatter(x, y, color="k")
plt.plot([0,x.max()], [0,a * x.max()])

ああああああああああ、全然やりたいことと違う
plotする時に、0からスタートすると駄目なのか。。

平均(mean)は、np.mean()だから、こうか？

x = np.array([1,2,3,4,5])
y = np.array([1571,1630,1681,1603,1623])

m = np.mean(y)

def reg1dim(x, y):
	a = np.dot(x, y)/ (x ** 2).sum()
	return a 

a = reg1dim(x, y)

plt.scatter(x, y, color="k")
plt.plot([0,x.max()], [m,a * x.max()])

あれ、なんかちゃうなー、 [m,a * x.max()]のところがおかしい。。

最小二乗法はこういう図になるはずなんだけど。。

最尤推定

最尤推定法: ある自由パラメーターが与えられた時に、モデルが実データを予測する確率を最大化するような自由パラメーターを求める方法
→　観測結果が偶然偏ってしまった場合や、適用する確率分布仮定を誤ると、見当違いな推定結果が導かれる

よって、回数を増やさないと、信用度が上がらない
性善説に立っているようですね。

最尤推定とは、パラーメータθに関する尤度関数L(θ)を最大化するθを求める
θを求めるには、一階級微分するが、対数尤度関数log[e]L(θ)を取る

公式
d/dθ*log[e]L(θ) = 0

最尤推定は、過去データから将来を予想する時に用いられることが多い
まぁ、データ量が上がれば精度も上がるということでしょう。

相関係数

nikkeiとdowの相関を見て見ましょう。
—
nikkei
2019/10/4 21,316 21,410 21,276 21,410
2019/10/3 21,422 21,437 21,277 21,341
2019/10/2 21,744 21,795 21,725 21,778
2019/10/1 21,831 21,938 21,811 21,885
2019/9/30 21,793 21,811 21,666 21,755
—
dow
2019年10月04日 26,573.72 26,271.70 26,590.74 26,271.70 224.49M 1.42%
2019年10月03日 26,201.04 26,039.02 26,205.20 25,743.46 249.02M 0.47%
2019年10月02日 26,078.62 26,425.86 26,438.04 25,974.12 312.73M -1.86%
2019年10月01日 26,573.04 26,962.54 27,046.21 26,562.22 272.18M -1.28%
2019年09月30日 26,916.83 26,852.33 26,998.86 26,852.33 228.48M 0.36%
—

import numpy as np
import pandas as pd

nikkei = [21410,21341,21778,21885,21755]
dow = [26573.72,26201.04,26078.62,26573.04,26916.83]
s1 = pd.Series(nikkei)
s2 = pd.Series(dow)
print(s1.corr(s2))

[vagrant@localhost python]$ python app.py
0.24461988452298797

あれ、思ったより全然相関してない??
dowに反応して日経が動く傾向とすれば、nikkeiを1日ずらせばいいのかな。

平均、分散、標準偏差

普通にやっても詰まらないので、日経平均の平均を取りましょう

日付、始値、高値、安値、終値

2019/10/4	21,316	21,410	21,276	21,410
2019/10/3	21,422	21,437	21,277	21,341
2019/10/2	21,744	21,795	21,725	21,778
2019/10/1	21,831	21,938	21,811	21,885
2019/9/30	21,793	21,811	21,666	21,755

nikkei = [21410,21341,21778,21885,21755]
print(sum(nikkei)/len(nikkei))

[vagrant@localhost python]$ python app.py
21633.8

分散は

import numpy as np

nikkei = [21410,21341,21778,21885,21755]
print(np.var(nikkei))

[vagrant@localhost python]$ python app.py
46880.56

標準偏差は、print(np.std(nikkei))とします。

[vagrant@localhost python]$ python app.py
216.5191908353622

あああ、今週は動きましたねーーーー

ロジスティック分布と正規分布

ロジスティック分布：連続確率分布の一つで、その累積分布関数がロジスティック関数であるもの？
パラメーターmu:平均
sigma:スケールパラメータ

import numpy as np 
import matplotlib.pyplot as plt 

sigma = 1.0
mu = 0

x = np.arange(-5., 5., 0.001)
y = np.exp(-(x-mu)/sigma) / (sigma*(1+np.exp(-(x-mu)/sigma))**2)
plt.plot(x,y)


正規分布と見た目は全く一緒だけど、ロジスティック分布は正規分布より裾野が長く、異なるようです。
– 密度関数は平均から離れても下がりにくい

あれ、でも正規分布の特徴である平均値、最頻値、中央値は一致してるのかな。

パレート分布

パレート分布とは所得者の分布？？

import numpy as np 
import matplotlib.pyplot as plt 
a, m = 6., 2.
s = (np.random.pareto(a, 1000) + 1) * m
count, bins, _ = plt.hist(s, 100, density=True)

fit = a * m ** a / bins ** (a + 1)
plt.plot(bins, max(count) * fit / max(fit), linewidth=2, color='r')

差も然りか。。

指数分布とは？

事象が連続して独立に一定の発生率で起こる過程
-> ランダムなイベントの発生間隔を表す分布
　e.g. 地震が起きる間隔、電球の寿命

f(x) = 1/μ e ^-x/μ (x>0)

import numpy as np 
import math
import matplotlib.pyplot as plt 

x = np.arange(-3,3.1,0.1)
y_2 = 2**x
y_3 = 3**x
y_e = math.e**x
y_0_5 = 0.5**x

plt.plot(x,y_2,label="a=2")
plt.plot(x,y_3,label="a=3")
plt.plot(x,y_e,label="a=e")
plt.plot(x,y_0_5,label="a=0.5")
plt.legend()
plt.savefig("01")

正規分布

正規分布:平均値の周辺にデータが集積するデータ分布
平均が50、対象が1000、標準偏差が10とする
確率密度を狭めていく

import numpy as np 
import matplotlib.pyplot as plt 

x = np.random.normal(50, 10, 1000)

plt.hist(x, bins=20)
plt.savefig("01")

plt.hist(x, bins=50)
plt.savefig("02")

plt.hist(x, bins=100)
plt.savefig("03")

まあ、正常なデータでしょうね。