import java.awt.RenderingHints;
import java.awt.geom.*;
import java.awt.Color;
import java.awt.BasicStroke;
public class tax31 extends JPanel{
public static void main(String[] args){
JFrame frame = new JFrame();
tax31 app = new tax31();
frame.getContentPane().add(app);
frame.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
frame.setBounds(0, 0, 500, 500);
frame.setTitle("Mirrlees(1971)の追試");
frame.setVisible(true);
}
public void paintComponent(Graphics g){
int n;
double data[]=new double[201];
Graphics2D g2 = (Graphics2D)g;
g2.setRenderingHint(RenderingHints.KEY_ANTIALIASING,
RenderingHints.VALUE_ANTIALIAS_ON);
data=tax();
for (n=1;n<99;n++){
g2.draw(new Line2D.Double(500*data[n+100], 500-500*data[n], 500*data[n+101],500-500*data[n+1]));
}
}
static double[] tax(){
double th[]= new double[101];
double c[]= new double[101];
double y[]= new double[101];
int opc[]= new int[101];
int opy[]= new int[101];
int opq[]= new int[101];
double u[][][]=new double[101][3][3];
double ww[][][]=new double[100][3][3];
double v[][][][]=new double[100][3][3][21];
int gotoc[][][][]=new int[100][3][3][21];
int gotoy[][][][]=new int[100][3][3][21];
int gotoq[][][][]=new int[100][3][3][21];
double endv[][]=new double[3][3];
int endc[][]=new int[3][3];
int endy[][]=new int[3][3];
int endq[][]=new int[3][3];
int s;
double ls;
double w;
double cs;
double ys;
double b1;
double b2;
double wel;
double maxwel;
int t;
int n;
int n1;
int n2;
int pp;
int q;
int j;
double c1;
double l1;
double h;
double u1;
int ns1;
int ns2;
int qx;
int nx1;
int nx2;
double vs;
double v1;
double maxv;
double dc;
double dy;
int e;
double data1[]=new double[101];
double data2[]=new double[201];
double data3[]=new double[201];
double maxy;
data1=ex1();
data2=ex2();
for (s=1;s<101;s++){
c[s]=data2[s];
y[s]=data2[s+100];
}
for (s=1;s<101;s++){
th[s]=data1[s];
}
h=0.001;
t=0;
while(t<1000){
for (s=1;s<101;s++){
for (n1=-1;n1<2;n1++){
for (n2=-1;n2<2;n2++){
c1=c[s]+n1*h;
l1=(y[s]+n2*h)/th[s];
u1=ux(c1,l1);
u[s][n1+1][n2+1]=u1;
}
}
}
for (s=1;s<100;s++){
for (n1=-1;n1<2;n1++){
for (n2=-1;n2<2;n2++){
c1=c[s]+n1*h;
l1=(y[s]+n2*h)/th[s+1];
u1=ux(c1,l1);
ww[s][n1+1][n2+1]=u1;
}
}
}
for (n1=-1;n1<2;n1++){
for (n2=-1;n2<2;n2++){
for (q=-10;q<11;q++){
v[1][n1+1][n2+1][q+10]=-999;
}
}
}
for (n1=-1;n1<2;n1++){
for (n2=-1;n2<2;n2++){
q=n2-n1;
v[1][n1+1][n2+1][q+10]=u[1][n1+1][n2+1];
}
}
for (s=2;s<100;s++){
for (n1=-1;n1<2;n1++){
for (n2=-1;n2<2;n2++){
for (q=-10;q<11;q++){
u1=u[s][n1+1][n2+1];
qx=q-n1+n2;
pp=0;
if (qx>10)pp=100;
if (qx<-10)pp=100;
if (pp>50)qx=0;
vs=-999;
ns1=0;
ns2=0;
for (nx1=-1;nx1<2;nx1++){
for (nx2=-1;nx2<2;nx2++){
v1=u1+v[s-1][nx1+1][nx2+1][qx+10];
if (ww[s-1][nx1+1][nx2+1]>u1)v1=-999;
if (v1>vs)ns1=nx1;
if (v1>vs)ns2=nx2;
if (v1>vs)vs=v1;
}
}
if (pp>50)vs=-999;
gotoc[s][n1+1][n2+1][q+10]=ns1;
gotoy[s][n1+1][n2+1][q+10]=ns2;
gotoq[s][n1+1][n2+1][q+10]=qx;
v[s][n1+1][n2+1][q+10]=vs;
}
}
}
}
for (n1=-1;n1<2;n1++){
for (n2=-1;n2<2;n2++){
u1=u[100][n1+1][n2+1];
qx=n2-n1;
vs=-999;
ns1=0;
ns2=0;
for (nx1=-1;nx1<2;nx1++){
for (nx2=-1;nx2<2;nx2++){
v1=u1+v[99][nx1+1][nx2+1][qx+10];
if (ww[99][nx1+1][nx2+1]>u1){
v1=-999;
}
if (v1>vs){
ns1=nx1;
}
if (v1>vs){
ns2=nx2;
}
if (v1>vs){
vs=v1;
}
}
}
endc[n1+1][n2+1]=ns1;
endy[n1+1][n2+1]=ns2;
endq[n1+1][n2+1]=qx;
endv[n1+1][n2+1]=vs;
}
}
maxv=-999;
ns1=0;
ns2=0;
for (n1=-1;n1<2;n1++){
for (n2=-1;n2<2;n2++){
if (endv[n1+1][n2+1]>maxv){
ns1=n1;
}
if (endv[n1+1][n2+1]>maxv){
ns2=n2;
}
if (endv[n1+1][n2+1]>maxv){
maxv=endv[n1+1][n2+1];
}
}
}
opc[100]=ns1;
opy[100]=ns2;
opc[99]=endc[opc[100]+1][opy[100]+1];
opy[99]=endy[opc[100]+1][opy[100]+1];
opq[99]=endq[opc[100]+1][opy[100]+1];
for (j=1;j<99;j++){
s=99-j;
opc[s]=gotoc[s+1][opc[s+1]+1][opy[s+1]+1][opq[s+1]+10];
opy[s]=gotoy[s+1][opc[s+1]+1][opy[s+1]+1][opq[s+1]+10];
opq[s]=gotoq[s+1][opc[s+1]+1][opy[s+1]+1][opq[s+1]+10];
}
e=0;
for (s=1;s<101;s++){
e=e+opc[s]*opc[s]+opy[s]*opy[s];
}
for (s=1;s<101;s++){
c[s]=c[s]+opc[s]*h;
y[s]=y[s]+opy[s]*h;
}
System.out.println(e);
if (e<2){
h=h/2;
}
if (h<0.00001){
t=10000;
}
t=t+1;
}
maxy=-999;
for (s=1;s<101;s++){
if (y[s]>maxy)maxy=y[s];
}
for (s=1;s<99;s++){
dc=c[s+1]-c[s];
dy=y[s+1]-y[s];
data3[s]=0;
if (dy>0)data3[s]=1-dc/dy;
}
for (s=1;s<100;s++){
data3[s+100]=y[s]/maxy;
}
return data3;
}
static double[] ex2(){
double th[]= new double[101];
int s;
double tl;
double tr;
double ls;
double w;
double cs;
double ys;
double b1;
double b2;
double tr1;
double tr2;
double tr3;
double wel;
double maxwel;
double maxtl;
double maxtr;
int t;
int n;
double h;
double c[]= new double[101];
double y[]= new double[101];
double data1[]=new double[101];
double data[]=new double[201];
int n1;
int n2;
data1=ex1();
for (s=1;s<101;s++){
th[s]=data1[s];
}
maxwel=-999;
maxtr=0;
maxtl=0;
for (n=20;n<45;n++){
tl=0.01*n;
tr1=0.01;
tr2=0.02;
tr=tr1;
b1=bud(tl,tr1,th);
t=0;
while (t<100) {
b2=bud(tl,tr2,th);
tr3=tr2-b2*(tr2-tr1)/(b2-b1);
tr1=tr2;
tr2=tr3;
b1=b2;
if (b2*b2<0.00001)t=1000;
t=t+1;
}
tr=tr2;
wel=seekwel(tl,tr,th);
if (wel>maxwel)maxtl=tl;
if (wel>maxwel)maxtr=tr;
if (wel>maxwel)maxwel=wel;
}
tl=maxtl;
tr=maxtr;
for (s=1;s<101;s++){
w=(1-tl)*th[s];
ls=(w-tr)/(2*w);
if (ls<0)ls=0;
c[s]=w*ls+tr;
y[s]=th[s]*ls;
}
for (s=1;s<101;s++){
data[s]=c[s];
data[s+100]=y[s];
}
return data;
}
static double ux(double c1,double l1){
double u1;
double c3;
double l3;
int pp;
pp=0;
c3=c1;
l3=l1;
if (c3<0)pp=100;
if (l3<0)pp=100;
if (l3>1)pp=100;
if (pp>50)c3=0.5;
if (pp>50)l3=0.5;
u1=Math.log(c3)+Math.log(1-l3);
if (pp>50)u1=-999;
return u1;
}
static double bud(double tl,double tr,double th[]){
double bx;
int s;
double w;
double y1;
double c1;
double l1;
bx=0;
for (s=1;s<101;s++){
w=(1-tl)*th[s];
l1=(w-tr)/(2*w);
if (l1<0)l1=0;
y1=th[s]*l1;
c1=w*l1+tr;
bx=bx+y1-c1;
}
return bx;
}
static double seekwel(double tl,double tr,double th[]){
double sw;
int s;
double w;
double l1;
double c1;
sw=0;
for (s=1;s<101;s++){
w=(1-tl)*th[s];
l1=(w-tr)/(2*w);
if (l1<0)l1=0;
c1=w*l1+tr;
sw=sw+ux(c1,l1);
}
return sw;
}
static double[] ex1(){
double p;
double mu;
double sig;
double yy;
double th[]= new double[101];
int s;
mu=0;
sig=0.39;
for (s=1;s<101;s++){
p=0.01*s-0.005;
yy=seeky(p,mu,sig);
th[s]=Math.exp(yy);
}
return th;
}
static double f(double x,double mu,double sig){
double pi,x1,x2,x3,fx;
pi = 3.1415;
x1=-Math.pow(x - mu,2) / (2*Math.pow(sig,2));
x2=Math.exp(x1);
x3=sig*Math.pow(2*pi,0.5);
fx=x2/x3;
return fx;
}
static double g(double y,double mu, double sig){
double gx,h,x;
int n,t;
gx=0;
h=0.001;
t=(int)(y/h);
for (n=-2000;n<t;n++){
x=n*h;
gx=gx+f(x,mu,sig)*h;
}
return gx;
}
static double seeky(double p,double mu,double sig){
double g1,g2,y1,y2,y3;
int t;
y1=0.4;
y2=-0.2;
g1=g(y1,mu,sig);
t=0;
while(t<100){
g2=g(y2,mu,sig);
y3=y2+(p-g2)*(y2-y1)/(g2-g1);
y1=y2;
y2=y3;
g1=g2;
if (Math.pow(p-g2,2)<0.0001)t=1000;
t=t+1;
}
return y2;
}
}
最終更新：2010年01月04日 17:00
