気象統計学 地人書館　鈴木栄一先生

#気象統計学　鈴木栄一先生　地人書館

peason_dis=function(y,a,b,alpha,beta){

val=log(gamma(a+b))-log(gamma(a))-log(gamma(b))+(a-1)*log(y-alpha)+(b-1)*log(beta-y)-(a+b+1)*log(beta-alpha)

return(exp(val))

}

x=rpois(100,8)

alphas=min(x);betas=max(x)

Z=(x-alphas)/(betas-alphas)

#pre-Broyden Algorithm p.211

f=function(s){

f1=s[1]/(s[1]+s[2])-mean(Z)

f2=s[1]*s[2]/((s[1]+s[2]+1)*(s[1]+s[2])^2)-var(Z)

return(c(f1,f2))

}

X=rep(1,2)

ite=100;h=0.01;eta=0.01

B=diag((f(X+h)-f(X))/h)

for(l in 1:ite){

if(sum(abs(f(X)))>10^(-9)){  

X_pre=X

X=X-solve(B)%*%f(X)

s=X-X_pre

q=f(X)-f(X_pre)

B=B+((q-B%*%s)/as.numeric(t(s)%*%s))%*%t(s)

print(f(X))

}}

plot(seq(min(x),max(x),0.001),peason_dis(seq(min(x),max(x),0.001),X[1],X[2],alphas,betas))

peason_dis_data=data.frame(x=seq(min(x),max(x),0.001),dis=peason_dis(seq(min(x),max(x),0.001),X[1],X[2],alphas,betas))




#超ガンマ分布

x=rpois(60,8);N=length(x);

h=0.01

X=rep(2,3)

eta=10^(-7)

ite=1000

log_G=sum(log(x))/N

for(l in 1:ite){

X_pre=X

alpha=X[1];beta=X[2];v=X[3]

mat=array(0,dim=c(3,3))

mat[1,1]=log_G+(v/(alpha^2))*(digamma(v/alpha+h)-digamma(v/alpha))/h

mat[1,2]=1/beta

mat[1,3]=-(digamma(v/alpha+h)-digamma(v/alpha))/(alpha*h)

mat[2,1]=beta*sum(x^alpha)+alpha*beta*sum(log(x)*x^alpha)

mat[2,2]=alpha*sum(x^alpha)

mat[2,3]=-N

mat[3,1]=beta*sum(log(x)*x^alpha)+alpha*beta*sum((log(x)^2)*x^alpha)

mat[3,2]=sum(log(x)*x^alpha)*alpha

mat[3,3]=-N*log_G

vec=c(-alpha*log_G+digamma(v/alpha)-log(beta),N*v-alpha*sum(x^alpha),N+N*v*log_G-alpha*beta*sum(log(x)*x^alpha))

vecs=mat%*%vec

X=X-eta*vecs

loglik=log(beta)*N*v/alpha+N*log(alpha)-N*log(gamma(v/alpha))-beta*sum(x^alpha)+(v-1)*sum(log(x))

print(loglik)

#print(sum(abs(X-X_pre)))

}







#ポリア・エッゲンベルガーの分布と台風襲来数の頻度に対する適合度検定

#台風襲来数
num=c(0:6)

#度数（年数）
freq=c(2,5,12,9,12,5,2)

#平均襲来数
ave=sum(num*freq/sum(freq))

#襲来数の分散
var=sum((num^2)*freq/sum(freq))-sum(num*freq/sum(freq))^2

h=ave;d=var/ave-1

polya_dis=function(s){

return(gamma(h/d+s)*((d/(1+d))^s)/(factorial(s)*gamma(h/d)*(1+d)^(h/d)))

}

pre=polya_dis(num)

fs=sum(freq)*pre

#ポリア・エッゲンベルガー分布はこのデータに適合しているか

kai2=sum((freq-fs)^2/fs)

qchisq(1-0.01,length(num)-2-1)

#p.88 母出現率に関する推測

#梅雨期の天気予報でAが441回中307回的中、Bが172回中94回的中した場合、どちらが的中させているか

p_A=307/441;p_B=94/172

k_A=307;N_A=441;k_B=94;N_B=172

alpha=0.25

p_A_upper=2*(k_A+1)*qf(1-alpha,2*(k_A+1),2*(N_A-k_A))/(2*(N_A-k_A)+2*(k_A+1)*qf(1-alpha,2*(N_A-k_A),2*(k_A+1)))

p_A_under=2*k_A/(2*k_A+2*(N_A-k_A+1)*qf(1-alpha,2*(N_A-k_A+1),2*k_A))

p_B_upper=2*(k_B+1)*qf(1-alpha,2*(k_B+1),2*(N_B-k_B))/(2*(N_B-k_B)+2*(k_B+1)*qf(1-alpha,2*(N_B-k_B),2*(k_B+1)))

p_B_under=2*k_B/(2*k_B+2*(N_B-k_B+1)*qf(1-alpha,2*(N_B-k_B+1),2*k_B))

#p_A_under>p_B_under?

static=(N_B-k_B)*k_A/((N_A-k_A+1)*(k_B+1))

qf(1-alpha,2*(k_B+1),2*(N_B-k_B))*qf(1-alpha,2*(N_A-k_A+1),2*k_A)


#p.224 censored sampling theoryと総雨量、最高風速

#A.K.Gupta

k=30;n=50

x=rpois(k,8)

d=max(x)-mean(x)

pthi=var(x)/(var(x)+d^2);p=k/n

#Z=1としてやってみる

Z=1

sigma_hat=d/Z;mu_hat=mean(x)+(sigma_hat^2-var(x))/d

log_lik=function(mu,sigma){

#val=sum(log(c(1:n)))-sum(log(c(1:(k-1))))-sum(log(c(1:(n-k))))-n*log(sigma/sqrt(2*pi))-sum((x-mu)^2)/(2*sigma^2)

val=-n*log(sigma/sqrt(2*pi))-sum((x-mu)^2)/(2*sigma^2)

fun=function(z){

return(exp(-(z-mu)^2/(2*sigma^2)))

}

val=val+(n-k)*log(integrate(fun,max(x),Inf)$value)

return((val))

}

h=0.01

v11=-(sigma_hat^2/n)*(log_lik(mu_hat+h,sigma_hat)-2*log_lik(mu_hat,sigma_hat)+log_lik(mu_hat-h,sigma_hat))/(h^2)

v22=-(sigma_hat^2/n)*(log_lik(mu_hat,sigma_hat+h)-2*log_lik(mu_hat,sigma_hat)+log_lik(mu_hat,sigma_hat-h))/(h^2)

v12=-(sigma_hat^2/n)*(log_lik(mu_hat+h,sigma_hat+h)-log_lik(mu_hat+h,sigma_hat)-(log_lik(mu_hat,sigma_hat+h)-log_lik(mu_hat,sigma_hat)))/(h^2)

sig_ij=(sigma_hat^2/n)*1/c(v11,v12,v22)


#鈴木栄一先生の方法

x=rpois(k,8)

sita_hat=k*((n-k)*max(x)+sum(x))^(-1)

V_sita_hat=sita_hat*(1-exp(-sita_hat*max(x)))^(-1)/n

E_x_n=mean(x)+(n-k)*max(x)/n

vec=c();count=0

while(count<(n-k))({

val=rexp(1,sita_hat)

if(val>max(x)){

count=count+1

vec=c(vec,val)

}

})

x_n_ave=mean(c(x,vec))

V_x_n=(1-exp(-mean(x)/x_n_ave))/(n*x_n_ave^2)



#予報の費用と的中率による最適化 p.263

N=10