へのフィット(続きの続きの続き)
そこで登場した、尤度、累積尤度をシミュレートして雰囲気をつかむためのエクセルはこちら。
デフォルトでは、1000個（N)の独立一様乱数(0-1)の最小値を100回サンプリング(計1000x100 乱数)し、その100個の最小値のセットから、もっともありそうなNの値を推定する(最尤推定量)とともに、尤度をN=nの関数として表示し、その尤度関数の累積分布をプロットするものである。
ファイルはこちら
Java ソース
package mlMinP;public class MLminP {
public double Nml;
public int k;
public double lnP;
public double[] as;
public double[] lnAs;
public double like0;
public double cumulLike0;
public double cumulLikeRatio0;
public double Lml;
public double lnLml;
/**
* @param args
*/
public static void main(String[] args) {
// TODO 自動生成されたメソッド・スタブ
int k=100;
int N=1000;
int dup=3;
double r=0.9;
double[] ps = new double[k];
for(int i=0;i<ps.length;i++){
//ps[i]=outMinP3(N,dup,r);
ps[i]=outMinP2(N,dup);
}
MLminP mlp = new MLminP(ps);
//mlp.PrintBriefMLminP("\t", "\n");
double ciIn = 0.95;
int type = 0;
double ceil = N*10;
double thres = 0.0000001;
CImlMinP ciM = new CImlMinP(mlp,type,ciIn,ceil,thres);
ciM.PrintCImlMinP("\t", "\n");
int type2 = 1;
CImlMinP ciM2 = new CImlMinP(mlp,type2,ciIn,ceil,thres);
ciM2.PrintCImlMinP("\t", "\n");
int numplot=100;
ciM.DrawPlot(numplot,"\t","\n");
double[] Nmlps = ciM.ObsQuantile(mlp,ps,mlp.Nml);
double[] ciLowps = ciM.ObsQuantile(mlp,ps,ciM.cilow);
double[] ciHighps = ciM.ObsQuantile(mlp,ps,ciM.cihigh);
for(int i=0;i<ps.length;i++){
System.out.println(ps[i]+"\t"+Nmlps[i]+"\t"+
ciLowps[i]+"\t"+ciHighps[i]);
}
StatUtils.MiscUtil.quickSortDouble(Nmlps, 0, Nmlps.length-1);
StatUtils.MiscUtil.quickSortDouble(ps, 0, ps.length-1);
for(int i=0;i<ps.length;i++){
double tmp=(double)(i+1)/(ps.length+1);
System.out.println(ps[i]+"\t"+Nmlps[i]+"\t"+
tmp);
}
}
public MLminP(double[] ps){
k=ps.length;
lnP=0;
for(int i=0;i<ps.length;i++){
lnP+=Math.log(1-ps[i]);
}
Nml=-(double)k/lnP;
as = new double[k+1];
lnAs = new double[k+1];
lnAs[k]=-Math.log(-lnP);
//as[k]=1/lnP;
for(int i=k-1;i>=0;i--){
lnAs[i]=lnAs[i+1]+Math.log(i+1)-Math.log(-lnP);
//as[i]=-Math.exp(lnAs[i]);
}
like0=0;
lnLml=k*Math.log(Nml)+(Nml-1)*lnP;
//Lml=Math.exp(lnLml);
//cumulLike0=-Math.exp(lnAs[0]-lnP);
cumulLikeRatio0=-Math.exp(lnAs[0]-1*lnP-lnLml);
//System.out.println("ratio0="+cumulLikeRatio0);
}
public double QuantileBothSide(MLminP mlp,double a,double ceil,double thres){
double ret = 0;
double equiv = mlp.seekUppEquiv(mlp, a,ceil, thres);
//System.out.println("equiv="+equiv);
double tmp =1-mlp.CumulStandLike(mlp, equiv);
//System.out.println("1="+mlp.CumulStandLike(mlp,a) + " 2="+ tmp);
ret = mlp.CumulStandLike(mlp,a) + 1-mlp.CumulStandLike(mlp, equiv);
return ret;
}
public double seekQuantileMLminPTwoSide(MLminP mlp,double q1, double ceil,double thres){
double ret = 0;
double q=1-q1;
/*
* 上側と下側を併せて1-q
* 上側の値に対する尤度と
* 下側の値に対する尤度とが同一であること
*/
double lower = 0;
//double higher = ceil;
double higher = Nml;//Nmlより下に下の閾値はあるので
double middle = higher/2;
double h=QuantileBothSide(mlp,higher,ceil,thres);
double l=QuantileBothSide(mlp,lower,ceil,thres);
double m=QuantileBothSide(mlp,middle,ceil,thres);
while(higher-lower>thres){
if(q>=l && q<=m){
higher=middle;
h=m;
middle=(lower+higher)/2;
m=QuantileBothSide(mlp,middle,ceil,thres);
}else if(q>=m && q<=h){
lower=middle;
l=m;
middle=(lower+higher)/2;
m=QuantileBothSide(mlp,middle,ceil,thres);
}else if(q>h){
higher*=2;
h=CumulStandLike(mlp,higher);
middle=(lower+higher)/2;
m=QuantileBothSide(mlp,middle,ceil,thres);
}
}
ret = (higher+lower)/2;
return ret;
}
/*
* Nmlより小さいnについて、そのnと同じ尤度を持つn'>Nmlを探す
*/
public double seekUppEquiv(MLminP mlp,double a, double ceil,double thres){
double ret = 0;
double lower = mlp.Nml;
double q = mlp.LikelihoodRatio(a);
//double higher = ceil;
double higher = ceil;//Nmlより下に下の閾値はあるので
double middle = (lower+higher)/2;
double l=mlp.LikelihoodRatio(higher);
double h=mlp.LikelihoodRatio(lower);
double m=mlp.LikelihoodRatio(middle);
while(higher-lower>thres){
if(q>=l && q<=m){
lower=middle;
h=m;
middle=(lower+higher)/2;
m=LikelihoodRatio(middle);
}else if(q>=m && q<=h){
higher=middle;
l=m;
middle=(lower+higher)/2;
m=mlp.LikelihoodRatio(middle);
}else if(q<l){
double tmphigh=higher*2;
//higher*=2;
lower=higher;
higher=tmphigh;
middle=(higher+lower)/2;
h=l;
m=mlp.LikelihoodRatio(middle);
l=mlp.LikelihoodRatio(higher);
}
}
ret = (higher+lower)/2;
return ret;
}
public double seekQuantileMLminP(MLminP mlp,double q,double ceil,double thres){
double ret = 0;
//double ceil = 10000000;
//double thres = 0.0000001;
double lower = 0;
double higher = ceil;
double middle = higher/2;
double h=CumulStandLike(mlp,higher);
double l=CumulStandLike(mlp,lower);
double m=CumulStandLike(mlp,middle);
while(higher-lower>thres){
if(q>=l && q<=m){
higher=middle;
h=m;
middle=(lower+higher)/2;
m=CumulStandLike(mlp,middle);
}else if(q>=m && q<=h){
lower=middle;
l=m;
middle=(lower+higher)/2;
m=CumulStandLike(mlp,middle);
}else if(q>h){
higher*=2;
h=CumulStandLike(mlp,higher);
middle=(lower+higher)/2;
m=CumulStandLike(mlp,middle);
}
}
ret = (higher+lower)/2;
return ret;
}
public double CumulativeLikeRatio(double n){
/*
* k*LN(N_ML)
* (N_ML-1)*LN(P)*/
double ret = 0;
if(n==0){
//for(int i=1;i<lnAs.length;i++){
//ret += -Math.exp(lnAs[i]+(n-1)*lnP-lnLml);
//}
//ret += -Math.exp(lnAs[0]-lnP);
ret += -Math.exp(lnAs[0]+(n-1)*lnP-lnLml);
//ret=0;
}else{
for(int i=0;i<lnAs.length;i++){
ret += -Math.exp(lnAs[i]+i*(Math.log(n))+(n-1)*lnP-lnLml);
}
}
return ret;
}
public double CumulativeLikeRatio0(double lnP, double lnLml){
/*
* k*LN(N_ML)
* (N_ML-1)*LN(P)*/
double ret = 0;
ret += -Math.exp(lnAs[0]-1*lnP-lnLml);
//ret=0;
return ret;
}
public double CumulStandLike(MLminP mlp,double n){
double ret = (mlp.CumulativeLikeRatio(n)-mlp.cumulLikeRatio0)/(-mlp.cumulLikeRatio0);
return ret;
}
public double CumulativeLike(double n){
double ret = 0;
if(n==0){
ret += -Math.exp(lnAs[0]-lnP);
}else{
for(int i=0;i<lnAs.length;i++){
ret += -Math.exp(lnAs[i]+i*Math.log(n)+(n-1)*lnP);
}
}
return ret;
}
public double CumulativeLike0(double lnP){
double ret = 0;
ret += -Math.exp(lnAs[0]-lnP);
return ret;
}
public double LikelihoodRatio(double n){
/*
* MLNに対する尤度比
*/
double ret = Likelihood(n)/Lml;
if(n==0){
ret =0;
}else{
ret = k*Math.log(n)+(n-1)*lnP-lnLml;
ret = Math.exp(ret);
}
return ret;
}
public double Likelihood(double n){
double ret = 0;
if(n==0){
}else{
ret = k*Math.log(n)+(n-1)*lnP;
ret = Math.exp(ret);
}
return ret;
}
public void PrintMLminP(String sep1,String sep2){
String ret = StringMLminP(sep1,sep2);
System.out.println(ret);
}public void PrintBriefMLminP(String sep1,String sep2){
String ret = StringBriefMLminP(sep1,sep2);
System.out.println(ret);
}
public String StringMLminP(String sep1,String sep2){
String ret ="";
ret += "Nml"+sep1+Nml + sep2;
ret += "k"+sep1+k + sep2;
ret += "lnP"+sep1+lnP + sep2;
//ret += "like0"+sep1+like0+sep2;
//ret += "cumulLike0"+sep1+cumulLike0+sep2;
//ret += "cumulLikeRatio0"+sep1+cumulLikeRatio0+sep2;
//ret += "LikeLM"+sep1+Lml+sep2;
//ret += "lnLikeLM"+sep1+lnLml+sep2;
ret += "lnAs"+sep2;
for(int i=0;i<lnAs.length;i++){
ret += lnAs[i]+sep1;
}
ret += sep2;
ret += "as"+sep2;
for(int i=0;i<as.length;i++){
ret += as[i]+sep1;
}
ret += sep2;
return ret;
}public String StringBriefMLminP(String sep1,String sep2){
String ret ="";
ret += "Nml"+sep1+Nml + sep2;
ret += "k"+sep1+k + sep2;
ret += "lnP"+sep1+lnP + sep2;
//ret += "like0"+sep1+like0+sep2;
//ret += "cumulLike0"+sep1+cumulLike0+sep2;
//ret += "cumulLikeRatio0"+sep1+cumulLikeRatio0+sep2;
//ret += "LikeLM"+sep1+Lml+sep2;
//ret += "lnLikeLM"+sep1+lnLml+sep2;
/*
ret += "lnAs"+sep2;
for(int i=0;i<lnAs.length;i++){
ret += lnAs[i]+sep1;
}
ret += sep2;
ret += "as"+sep2;
for(int i=0;i<as.length;i++){
ret += as[i]+sep1;
}
ret += sep2;
*/
return ret;
}
public static double outMinP(int N){
double ret = Math.random();
for(int i=1;i<N;i++){
double tmp = Math.random();
if(tmp<ret){
ret=tmp;
}
}
return ret;
}
public static double outMinP2(int N,int dup){
double ret = Math.random();
int count=1;
for(int i=1;i<N;i++){
double tmp=ret;
if(count==0){
tmp = Math.random();
}
if(count==dup){
count=0;
}else{
count++;
}
//double tmp = Math.random();
if(tmp<ret){
ret=tmp;
}
}
return ret;
}
public static double outMinP3(int N,int dup,double r){
double ret = Math.random();
double prev=ret;
int count=1;
double tmp=Math.random();;
int df=1;//自由度
for(int i=1;i<N;i++){
if(count==0){
tmp = Math.random();
}else{
try{
//double a = chidf1.cumulativeProbability(6.635);
org.apache.commons.math.distribution.NormalDistributionImpl chidf1 =
new org.apache.commons.math.distribution.NormalDistributionImpl(tmp,tmp*(1-r));
double b = chidf1.inverseCumulativeProbability(Math.random());
//System.out.println("a="+a);
//System.out.println("b=\t"+b);
if(b>0 && b<1){
tmp=b;
}
} catch (org.apache.commons.math.MathException ex) {
}
}
if(count==dup){
count=0;
}else{
count++;
}
if(tmp<ret){
ret=tmp;
}
}
return ret;
}
}
package mlMinP;public class CImlMinP {
public MLminP mlp;
public int cimethod;//0 both 1 even
public double cirange;
public double cilow;
public double cihigh;
public double[] cutoff = {0.01,0.05,0.1};
public double[][] nominalP = new double[3][cutoff.length];
public double ceil;
public double thres;
public CImlMinP(MLminP mlp_, int cimethod_,double cirange_,double ceil_,double thres_){
mlp=mlp_;
cimethod=cimethod_;
cirange=cirange_;
ceil=ceil_;
thres=thres_;
if(cimethod==0){
cilow=mlp.seekQuantileMLminPTwoSide(mlp,cirange,ceil,thres);
cihigh=mlp.seekUppEquiv(mlp, cilow, ceil, thres);
}else if(cimethod==1){
double low = (1-cirange)/2;
double high = 1-low;
cilow=mlp.seekQuantileMLminP(mlp,low,ceil,thres);
cihigh=mlp.seekQuantileMLminP(mlp,high,ceil,thres);
}
nominalP = new double[3][cutoff.length];
for(int i=0;i<nominalP[0].length;i++){
nominalP[0][i]=1-Math.pow((1-cutoff[i]),(1/mlp.Nml));
nominalP[1][i]=1-Math.pow((1-cutoff[i]),(1/cilow));
nominalP[2][i]=1-Math.pow((1-cutoff[i]),(1/cihigh));
}}
public double[] ObsQuantile(MLminP mlp,double[] ps,double n){
double[] ret = new double[ps.length];
for(int i=0;i<ret.length;i++){
ret[i]=1-Math.pow((1-ps[i]),n);
}
return ret;
}
public String StringCImlMinP(String sep1,String sep2){
String ret ="";
ret += mlp.StringBriefMLminP(sep1,sep2);
ret += sep2+"method"+sep1+cimethod+sep2;
ret += "cirange"+sep1+cirange+sep2;
ret += "Nml"+sep1+mlp.Nml+sep2;
ret += sep1+sep1+"like"+sep1+"cumul"+sep1;
for(int i=0;i<cutoff.length;i++){
ret += cutoff[i]+sep1;
}
ret += sep2;
for(int i=0;i<nominalP.length;i++){
if(i==0){
ret += "nML"+sep1+mlp.Nml+sep1;
double tmp = mlp.LikelihoodRatio(mlp.Nml);
ret += tmp+sep1;
tmp = mlp.CumulStandLike(mlp, mlp.Nml);
ret += tmp+sep1;
}else if(i==1){
ret += "ciLow"+sep1+cilow+sep1;
double tmp = mlp.LikelihoodRatio(cilow);
ret += tmp+sep1;
tmp = mlp.CumulStandLike(mlp,cilow);
ret += tmp+sep1;
}else if(i==2){
ret += "ciHigh"+sep1+cihigh+sep1;
double tmp = mlp.LikelihoodRatio(cihigh);
ret += tmp+sep1;
tmp = mlp.CumulStandLike(mlp,cihigh);
ret += tmp+sep1;
}
for(int j=0;j<nominalP[i].length;j++){
ret += nominalP[i][j]+sep1;
}
ret += sep2;
}
return ret;
}
public void PrintCImlMinP(String sep1,String sep2){
String ret=StringCImlMinP(sep1,sep2);
System.out.println(ret);
}
public void DrawPlot(int numplot,String sep1,String sep2){
double nlen=(double)numplot;
double ninit = cilow-(mlp.Nml-cilow)*2;
if(ninit<0){
ninit=0;
}
double nend = cihigh+(cihigh-mlp.Nml)*2;
double nkizami = (nend-ninit)/nlen;
for(int i=0;i<numplot;i++){
double tmpn=ninit+i*nkizami;
String tmp =tmpn+sep1;
//double like = mlp.Likelihood(tmpn);
double likeRatio = mlp.LikelihoodRatio(tmpn);
//System.out.println("like="+like);
//System.out.println("likeratio="+likeratio);
tmp+=likeRatio+sep1;
//double cumlike = mlp.CumulativeLike(tmpn);
//double cumulLikeRatio = mlp.CumulativeLikeRatio(tmpn);
double cumlikeRatioStand = mlp.CumulStandLike(mlp,tmpn);
//System.out.println("cumlike="+cumlike);
//System.out.println("cumlikeRatio="+cumlikeRatio);
tmp+=cumlikeRatioStand+sep1;
System.out.println(tmp);
}
/*
*
int nlen = 1000;
double nkizami = 400;
double ninit = 50000;
double[] ns = new double[nlen];
for(int i=0;i<ns.length;i++){
ns[i]=ninit+i*nkizami;
}
for(int i=0;i<ns.length;i++){
String tmp =ns[i]+"\t";
double like = mlp.Likelihood(ns[i]);
double likeRatio = mlp.LikelihoodRatio(ns[i]);
//System.out.println("like="+like);
//System.out.println("likeratio="+likeratio);
tmp+=like+"\t"+likeRatio+"\t";
double cumlike = mlp.CumulativeLike(ns[i]);
double cumulLikeRatio = mlp.CumulativeLikeRatio(ns[i]);
double cumlikeRatioStand = mlp.CumulStandLike(mlp,ns[i]);
//System.out.println("cumlike="+cumlike);
//System.out.println("cumlikeRatio="+cumlikeRatio);
tmp+=cumlike+"\t"+cumulLikeRatio + "\t"+cumlikeRatioStand+"\t";
System.out.println(tmp);
}
*/
}
}
