養花種魚數月亮賞星星

這是留言板-既然來了，就留個言再走吧！

2019-07-22T14:43:00.014+08:00

Please leave your name and tell me what you thought about this site. Comments, suggestions and views are welcomed.

Mathematica 教學：對遺失資料做線性插值

2012-02-26T15:37:00.002+08:00

datatemp={
{-818.2106,-6755.208,207.9009},{0,0,0},{0,0,0},{0,0,0},
{0,0,0},{0,0,0},{0,0,0},{0,0,0},{0,0,0},{0,0,0},{0,0,0},
{0,0,0},{0,0,0},{0,0,0},{0,0,0},{0,0,0},{0,0,0},{0,0,0},
{0,0,0},{0,0,0},{0,0,0},{0,0,0},{0,0,0},{0,0,0},{0,0,0},
{0,0,0},{0,0,0},{-767.4377,-6709.237,251.7202},
{-768.7812,-6700.438,253.0015},{-768.6588,-6701.238,252.8851},
{-767.7062,-6697.338,254.8504},{-765.7729,-6696.207,257.0257},
{-766.2335,-6689.349,259.0551},{-761.753,-6693.4,261.8711},
{-762.6193,-6685.736,264.3021},{-757.8909,-6690.449,267.0908},
{-754.3521,-6692.645,269.7541},{-752.003,-6692.321,272.2921},
{-750.8434,-6689.48,274.7048},{-750.8735,-6684.12,276.9922},
{-749.6801,-6681.277,279.5483},{-748.7109,-6677.929,282.1369},
{-747.9661,-6674.078,284.7577},{-747.4456,-6669.722,287.4109},
{-747.1495,-6664.861,290.0965},{-744.986,-6667.709,291.8721},
{0,0,0},{0,0,0},{0,0,0},{0,0,0},{0,0,0},{0,0,0},{0,0,0},
{0,0,0},{0,0,0},{0,0,0},{0,0,0},{0,0,0},{0,0,0},{0,0,0},
{0,0,0},{0,0,0},{0,0,0},{0,0,0},{0,0,0},{0,0,0},
{-714.5345,-6624.671,368.8935}
}

mytrans[datatemp_List]:=Block[{data,zeroindex},

(*先將讀入資料加進流水號*)

data=Flatten[{#,datatemp[[#]]}]&/@Range[Length@datatemp];

(*列出斷點*)

zeroindex=Split[Select[data,#[[2]]==0&][[All,1]],#2-#1==1&];

(*斷點內線性插值*)

test[z_]:=Block[{incr},

(*找出斷點前後的資料*)

data[[#]]&/@{zeroindex[[z,1]]-1,zeroindex[[z,-1]]+1};

(*計算插值並替換資料*)

incr=MapThread[#2-#1&,data[[#]]&/@{zeroindex[[z,1]]-1,
zeroindex[[z,-1]]+1}];
Map[(data[[#]]=data[[#-1]]+incr/incr[[1]])&,zeroindex[[z]]]];

(*計算所有斷點插值*)

test[#]&/@Range[Length@zeroindex];
data
]


With[{dd=mytrans[datatemp]},
GraphicsGrid@
Transpose@{ListPlot[datatemp[[All,#]],PlotRange->All]&/@{1,2,3},
ListPlot[dd[[All,#]],PlotRange->All]&/@{2,3,4}}]

Mathematica 教學：利用 Google Scholar 擷取期刊引用數

2012-02-22T11:55:00.005+08:00

(* 指定期刊在Google Scholar上的資訊量 *)

num[journal_,year_]:=ToExpression[StringReplace[Cases[Import["http://scholar.google.com/scholar?hl=en&num=100&q=&as_publication=%22"<>StringReplace[ToString[journal]," "->"+"]<>"+%22&btnG=Search&as_sdt=0%2C5&as_ylo="<>ToString[year]<>"&as_yhi="<>ToString[year]<>"&as_vis=0","XMLObject"],{"Results ",XMLElement["b",{},{"1"}],x___}:>x,Infinity][[4,3,1]],","->""]];

(* 建立超連結 *)

link[journal_,year_,i_]:="http://scholar.google.com/scholar?start="<>ToString[i]<>"hl=en&num=100&q=&as_publication=%22"<>StringReplace[ToString[journal]," "->"+"]<>"+%22&btnG=Search&as_sdt=0%2C5&as_ylo="<>ToString[year]<>"&as_yhi="<>ToString[year]<>"&as_vis=0"
temp[journal_String,year_]:=Block[{mynum,dd1,dd2,dd3,length,aa,tempdata1,tempdata2,t1,t2,rule,title,author,cite},
mynum=100(Floor[num[journal,year]/100]);

(* 待刪除的XML引數 *)

rule=XMLElement["span",___]|XMLElement["input",___]|XMLElement["td",___]|XMLElement["br",___]|XMLElement["table",___]|XMLElement["div",{"class"->"gs_rs"},___]|XMLElement["form",___]|{"shape"->"rect","href"->___};

(* 剔除超連結非來自期刊公司資料庫的指定期刊在Google Scholar上的引用數 *)

tempdata1=Cases[Flatten[Table[(
dd1=Import[link[journal,year,i],"XMLObject"];
length=Length[dd1[[2,3,2,3]]];
Table[
(dd2=If[Length@Flatten@StringPosition[Cases[dd1[[2,3,2,3]][[z]],XMLElement["a",{"shape"->"rect","href"->hyperlink_},__]:>hyperlink,Infinity],"sciencedirect"]==0,{},If[Length[Cases[dd1[[2,3,2,3]][[z]],XMLElement["span",{"class"->"gs_ctu"},{"[CITATION]"}]|XMLElement["span",{"class"->"gs_ctc"},{"[PDF]"}]|XMLElement["span",{"class"->"gs_ctc"},{"[HTML]"}],Infinity]]==1,{},dd1[[2,3,2,3]][[z]]]];
dd3=DeleteCases[dd2,rule,Infinity]/.XMLElement->List;
title=Cases[dd3,{_,{"class"->"gs_rt"},title_}:>title,Infinity];
author=Cases[dd3,{_,{"class"->"gs_a"},{author_,___}}:>author,Infinity];
cite=Cases[dd3,{_,{"class"->"gs_fl"},cite_}:>cite,Infinity];
Flatten[{title,author,cite},1]
),{z,length}]),{i,0,mynum,100}],1],{_,___}];
tempdata2=Table[{If[Length[tempdata1[[z]][[1]]]==1,tempdata1[[z]][[1,1,-1]],tempdata1[[z]][[1,1,-1]]<>" "<>tempdata1[[z]][[1,2,-1]]],tempdata1[[z]][[2]],If[Length[tempdata1[[z]]]<3,{},tempdata1[[z]][[3,1,2]]]},{z,Length@tempdata1}] ]

(* 範例：將European Journal of Operational Research在2008年的文章資訊匯出到test.xls *)

Export["test.xls", temp["European Journal of Operational Research", 2008]]  

(*將European Journal of Operational Research在2000~2010年的所有文章資訊輸出*)  

ejordata=Block[{data,cite}, data=temp["European Journal of Operational Research",#]; cite=Total@Table[ If[Length[Flatten[StringPosition[data[[z,-1]],"Cited"]]]>0,
Read[StringToStream@data[[z,-1,1]],{Word,Word,
Number}][[-1]],0],{z,Length@data}];
{Length@data,cite,cite/Length[data]}
]&/@Range[2000,2010]

(* EJOR各年度在 Google Scholar 的引用量 *)

BarChart[ejordata[[All, 2]], ChartLabels -> Range[2000, 2010, 1], ImageSize -> 500]

花海

2012-02-05T01:01:00.004+08:00

最近高雄市政府在美術館這邊設立了很多簡易公園，原本想說何謂簡易公園，原來是一片花海。人在花海中散步，也不自覺的翩翩起舞～

Mathematica 教學：檔案複製、搬移

2012-02-04T13:45:00.003+08:00

最近寫一本講義，用到的圖檔實在太多，整理資料實在費時。所以就利用Mathematica來做一些資料搬移複製的工作。

(*檔案位置*)
path = NotebookDirectory[];
(*設定目錄位置*)
SetDirectory[path];
(*輸入LaTeX檔Compile後的log檔*)
test = Import["Ch1.log"];
(*抓出該TeX檔所使用的所有eps檔*)
epsdata=StringTrim@Flatten[StringCases[#,x__~~".eps    Graphic file (type eps)"->x]&/@
Select[test,Length[#]==1&]]
(*新增一個資料夾CH1，將Ch1.tex所使用到的eps全部複製到Ch1*)
CreateDirectory["Ch1]
CopyFile[path<>#<>".eps", "path"<>"Ch1/"<>#<>".eps"]&/@epsdata

Mathematica 教學：World life data with mixed statistical graphics

2012-02-01T00:44:00.010+08:00

World life data with mixed statistical graphics

data=Import["https://www.cia.gov/library/publications/the-world-factbook/fields/2102.html","Data"];

datatemp=Select[data[[5;;-3]],Length[Tally[StringPosition[#,"NA"]]]==1&];

yearsget[cdata_]:=Block[{cc=cdata,strpos},
strpos=StringPosition[cc,{"total population:","years male:","years female:"}];
StringTake[cc,{{strpos[[1,2]]+1,strpos[[2,1]]-1},{strpos[[2,2]]+1,
strpos[[3,1]]-1},{strpos[[3,2]]+1,strpos[[3,2]]+5}}]
]

lifedata=Map[Flatten@{StringTake[#[[1]],{1,-2}],
ToExpression@yearsget[#[[2]]]}&,datatemp]

Needs["ComputationalGeometry`"]

SortBy[lifedata,Last]//TableForm

leq[ptstemp_List]:=Block[{pts=ptstemp,slope,x,y},
If[Abs[pts[[2,1]]-pts[[1,1]]]<10^-9,x==pts[[1,1]], y==(pts[[2,2]]-pts[[1,2]])/(pts[[2,1]]-pts[[1,1]])*(x- pts[[1,1]])+pts[[1,2]]]]; bagscheme[{sourcetemp_List,ptstemp_List}]:= Block[{source=sourcetemp,pts=ptstemp,ch,temp,newlineeq, newpts,x,y,newch,newchpts}, ch=Quiet@ConvexHull@Flatten[{source,pts},1]; temp=source[[#]]&/@Select[ch,#<=Length@source&]; source= source[[#]]&/@ Complement[Range[Length@source], Select[ch,#<=Length@source&]];  newlineeq=Flatten[{{{#[[1]],#[[-2]]},{#[[-1]],#[[2]]}},{{#[[1]],#[[ 3]]},{#[[-1]],#[[2]]}}}&/@ Table[RotateLeft[temp,i],{i,0,Length@temp-1}],1]; newpts={x,y}/.Solve[leq[#]&/@#,{x,y}][[1]]&/@newlineeq; N@{source,newpts}];  mybag=NestWhileList[bagscheme,{lifedata[[All,{3,4}]],{}}, Length@#[[1]]>Length[lifedata[[All,{3,4}]]]/2&,1];
mybagcontour[i_]:=
Tally[Flatten[mybag[[-i]],1]][[All,1]][[#]]&/@(Quiet@
ConvexHull@Tally[Flatten[mybag[[-i]],1]][[All,1]]);

Grid[{{Graphics[{Opacity[0.8],Green,
Cases[BoxWhiskerChart[lifedata[[All,4]],"Notched"],
GraphicsGroupBox[box___]:>box,Infinity]/.
PolygonBox[a__]:>PolygonBox[Evaluate@(If[#[[1]]>
1,{#[[1]]+1,#[[2]]},{#[[1]]-1,#[[2]]}]&/@
a)]},PlotRange->{Automatic,{20,100}},ImageSize->{50,700}],
ListPlot[lifedata[[All,{3,4}]],AspectRatio->1.5,
ImageSize->{500,700},
Axes->False,Frame->True,
GridLines->{Range[20,100,2],Range[20,100,2]},
AxesOrigin->{20,20},
PlotRange->{{20,100},{20,100}},
PlotStyle->Black,
Epilog->{
ConvexHullpoint=Quiet@ConvexHull[lifedata[[All,{3,4}]]];
Opacity[0.75],EdgeForm[None],
Red,
Cases[Histogram[lifedata[[All,3]],{25,100,5}],
RectangleBox[a__,b__,c___]:>
Rectangle[
a+{0,20},{b[[1]],20+a[[2]]+(b[[2]]-a[[2]])/2}],
Infinity],
Green,
Cases[
Histogram[lifedata[[All,4]],{25,100,5},BarOrigin->Left],
RectangleBox[a__,b__,c___]:>
Rectangle[
a+{20,0},{(b[[1]]-a[[1]])/2+a[[1]]+20,b[[2]]}],
Infinity],
Opacity[1],
Green,Thickness[0.01],
Line[{#,
LinearModelFit[lifedata[[All,{3,4}]],x,x][#]}&/@{20,
100}],
Brown,Thickness[0.01],
Line[{LinearModelFit[lifedata[[All,{4,3}]],y,
y][#],#}&/@{20,100}],
Black,
{Text[lifedata[[#,1]],lifedata[[#,{3,4}]]],
Point[lifedata[[#,{3,4}]]]}&/@ConvexHullpoint,
Blue,Opacity[0.25],EdgeForm[Blue],
Polygon[lifedata[[#,{3,4}]]&/@ConvexHullpoint],
Red,Opacity[0.25],EdgeForm[Red],
Polygon[Insert[mybagcontour[1],mybagcontour[1][[1]],-1]],Red,
Opacity[1],PointSize[0.0125],
Point@Quiet@ConvexHullMedian@lifedata[[All,{3,4}]]
}
]},
{"",Graphics[{Opacity[0.75],Red,
Cases[BoxWhiskerChart[lifedata[[All,3]],"Notched",
BarOrigin->Left],GraphicsGroupBox[box___]:>box,
Infinity]/.
PolygonBox[a__]:>
PolygonBox[
Evaluate@(If[#[[2]]>
1,{#[[1]],#[[2]]+1},{#[[1]],#[[2]]-1}]&/@
a)]},PlotRange->{{20,100},Automatic},
ImageSize->{500,50}]}}]

Mathematica 教學：表格自動取小數點及補零

2012-01-08T00:16:00.004+08:00

Mathematica可以將表格轉換成標準的LaTeX格式的表格。但是唯一的缺點是要自己條小數點，每次一個一個調整尤其是在每個欄位要求又不一樣時實在煩人。所以寫了以下這個函數，會根據指定的要求自動將小數點四捨五入或補零。懶得話，直接將程式加進init.m，每次程式打開會自動載入～～懶，的確是一種美德！

(*init.m路徑*)
$UserBaseDirectory <> "/Kernel/init.m"

data={{9549.26,9548.13,9549.12,0.351504,182.612,107.94},{9666.36,
9654.91,9665.92,2.17506,138.135,102.33},{8917.18,8917.18,
8917.18,1.45152*10^-7,75.7518,99.75},{8917.18,8917.18,8917.18,
1.22834*10^-10,340.693,472.25}};

mynumform[data_List,decimal_List]:=Block[{d=decimal,mydata=data},
Table[Map[NumberForm[#[[1]],{10,#[[2]]},
(*在-Infinity跟Infinity之間不用以指數形式表示*)
ExponentFunction->(If[-10^100<#<10^100,Null,#]&)]&,
Transpose@{mydata[[z]],d}],{z,Length@mydata}]
]

TableForm@mynumform[data, {2, 2, 2, 4, 3, 3}]

Mathematica 教學消除Piecewise函數在x軸上的零值

2011-12-16T19:02:00.007+08:00

利用Piecewise所建立的函數在x軸未定義的地方會顯示零值


f[x_]:=Piecewise[{{x^2,x > 1},{x,0< x <1},{-x,-1< x <0}}]

Plot[f[x],{x,-2,2},PlotStyle->Thickness[0.02]]

(*將y=0的點剔除*)
Plot[f[x],{x,-2,2},PlotStyle->Thickness[0.02]]/.Line[pt___]:>Line[Select[pt,#[[2]]!=0&]]

Mathematica 教學 ListPlot without PlotLegend

2011-12-14T20:21:00.006+08:00

寫了一個新版本，可以支援Filling這個Option。圖形的表現雖不典雅，但堪用！科科～

mylegend[plot_Graphics,legend_List]:=
Block[{p=plot,l=legend,colortemp,color,Opacitytemp,temp},
(*擷取色彩*)
colortemp=Cases[p,Hue[a_,b_,c_]:>Hue[a,b,c],Infinity];
color=colortemp[[-Length@l;;-1]];
(*配對圖例及色彩*)
temp={color[[#]],l[[#]]}&/@Range[Length@color];
(*建立圖例*)
Labeled[p,
Grid[{Graphics[{#[[1]],Thickness[0.1],Opacity[0.75],
Line[{{0,0},{1,0}}]},ImageSize->{24,24},
AspectRatio->8/24,ImagePadding->0],#[[2]]}&/@
temp],
(*指定圖例位置*)
{{Right,Top}}]]

(*Example*)
mylegend[ListLinePlot[Table[Accumulate[RandomReal[{-1,1},250]],{3}],
Filling->0],{"a","b","c"}]

mylegend[ListLinePlot[Table[Accumulate[RandomReal[{-1,1},250]],{10}]],
CharacterRange["a","z"][[1;;10]]]

上面的程式使用系統預設的色彩，下面程式可以自行指定顏色的話

mylegend[plot_Graphics,legend_List]:=
Block[{p=plot,l=legend,colortemp,color,Opacitytemp,temp},
(*擷取色彩*)
(*判斷是否自行者定顏色*)
colortemp=If[Cases[p,RGBColor[a_,b_,c_]:>RGBColor[a,b,c],
Infinity]=={},
Cases[p,Hue[a_,b_,c_]:>Hue[a,b,c],Infinity],
Cases[p,RGBColor[a_,b_,c_]:>RGBColor[a,b,c],Infinity]];
color=colortemp[[-Length@l;;-1]];

(*配對圖例及色彩*)
temp={color[[#]],l[[#]]}&/@Range[Length@color];

(*建立圖例*)Labeled[p,
Grid[{Graphics[{#[[1]],Thickness[0.1],Opacity[0.75],
Line[{{0,0},{1,0}}]},ImageSize->{24,24},
AspectRatio->8/24,ImagePadding->0],#[[2]]}&/@
temp],
(*指定圖例位置*){{Right,Top}}]
]
mylegend[ListLinePlot[Table[Accumulate[RandomReal[{-1,1},250]],{3}],
Filling->0,PlotStyle->{Red,Blue,Green}],
{"a","b","c"}]

Mathematica 教學 Plot 圖形旋轉

2011-12-12T21:48:00.002+08:00

g=Plot[10Sin[x]/x,{x,0,10}]

(*將圖形座標取出並旋轉，接著再指定顏色*)
f[scale_,theta_,color_]:=g/.{Line[a___]:>Line[scale*RotationMatrix[thetaDegree].#&/@a],
Hue[___]:>color}

Manipulate[
Show[g,f[s,t,Green],PlotRange->{{-20,20},{-20,20}},
AspectRatio->1],{{s,1,"Scale"},0,2},{{t,45,"Degree"},0,360}]

Mathematica 教學 Logistic regression

2011-12-10T23:56:00.005+08:00

Logistic regression

data={{1,0},{1,1},{1,1},{2,0},{2,0},{2,1},{2,1},{2,1},{2,1},{3,1},
{3,1},{3,1},{3,1},{3,1}};

(* 成功機率 *)
pi[x_]:=1/(1+Exp[-(a+b*x)])

(* Likelihood function*)
obj=Apply[Times,pi[#[[1]]]^#[[2]]*(1-pi[#[[1]]])^(1-#[[2]])&/@data]

(* 求解a and b *)

FindMaximum[Log[obj],{a,0.1},{b,0.1}]
Plot[pi[x]/.FindMaximum[obj,{a,b}][[2]],{x,1,3},PlotRange->All]

LaTeX titlesec 模版

2011-12-10T16:31:00.006+08:00

無聊修改一下講義使用的章節標題格式，使用titlesec這個Package。感覺這樣比較像一本書~~

\definecolor{titlecolor}{rgb}{0.503906,0.503906,0.734375}
\newcommand{\hwyk}{\CJKfamily{cwhbb}}

\newcommand{\mytitle}[1]{
\begin{tabular}{p{0.01\textwidth}p{0.99\textwidth}}
\cellcolor{black} &\cellcolor{titlecolor} \textcolor{white}
{\newline\hwyk \LARGE 第\ \thechapter  \ 章 \  \ #1}
\end{tabular}
\arrayrulewidth=0.4pt
}

\newcommand{\mysection}[1]{
\setlength\arrayrulewidth{1pt}\arrayrulecolor{titlecolor}
\begin{tabular}{p{0.01\textwidth}p{0.99\textwidth}}
\hline
\cellcolor{titlecolor} &  \textcolor{black}{\hwyk \LARGE \thesection ~ #1}
\end{tabular}
\arrayrulewidth=0.4pt
}

\definecolor{titlecolor}{rgb}{0.503906,0.503906,0.734375}
\newcommand{\hwyk}{\CJKfamily{cwhbb}}


\titleformat{\chapter}[hang]{\hwyk \LARGE   \sf}
{}{0mm}{\hspace{-0.4cm}\mytitle}

\titleformat{\section}[hang]{\hwyk \LARGE   \sf}
{}{0mm}{\hspace{-0.5cm}\mysection}

這種開場白真的有人願意幫忙填問卷嗎？

2011-12-09T19:57:00.002+08:00

這種開場白真的有人願意幫忙填問卷嗎？改成正妹付照片比較實在。

Mathematica 教學 Kernel Fisher discriminant analysis

2011-12-09T00:04:00.008+08:00

Kernel Fisher discriminant analysis

2-Dimension

data=SplitBy[Import["http://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data"],Last];

X=Flatten[data[[{1,3}]],1][[All,{1,2}]];
L=Length[X];
L1=50;L2=50;

(* Kernel function *)
RBFKernel[x_,y_,\[Gamma]_]:=Exp[-\[Gamma](x-y).(x-y)]

(*  Kernel matrix  *)
ttX=Table[RBFKernel[X[[i]],X[[j]],2],{i,100},{j,100}];

(* Sum square of Total *)
SST=Covariance[ttX]*(L-1);

(* Sum square of  within-class *)
SSW=(L1-1)*Covariance[ttX[[1;;50]]]+(L2-1)*Covariance[ttX[[51;;100]]];

(* Sum square of  between-class  *)
SSB=SST-SSW;

(* Eigen system *)
e1=Eigensystem[Inverse[SSW+IdentityMatrix[L]*L].SSB];

(* Kernel Fisher Discriminant function *)

kfd1=e1[[2,1]].Table[RBFKernel[X[[i]],{x1,x2},2],{i,L}];
P1=kfd1/.{x1->X[[#,1]],x2->X[[#,2]]}&/@Range[50];
P2=kfd1/.{x1->X[[#,1]],x2->X[[#,2]]}&/@Range[51,100,1];
kfd1=kfd1-(Mean[P1]StandardDeviation[P2]+
Mean[P2]StandardDeviation[P1])/(StandardDeviation[P2]+StandardDeviation[P1]);

(* Plot *)
Show[
ContourPlot[kfd1==0,{x1,4,8},{x2,2,5}],
ListPlot[{X[[1;;50]],X[[51;;100]]},PlotStyle->PointSize[0.02]]
]

3-Dimension

X=Flatten[data[[{1,3}]],1][[All,{1,2,3}]];
ttX=Table[RBFKernel[X[[i]],X[[j]],2],{i,100},{j,100}];

SST=Covariance[ttX]*(L-1);
SSW=(L1-1)*Covariance[ttX[[1;;50]]]+(L2-1)*Covariance[ttX[[51;;100]]];
SSB=SST-SSW;

e1=Eigensystem[Inverse[SSW+IdentityMatrix[L]*L].SSB];

kfd1=e1[[2,1]].Table[RBFKernel[X[[i]],{x1,x2,x3},2],{i,L}];
P1=kfd1/.{x1->X[[#,1]],x2->X[[#,2]],x3->X[[#,3]]}&/@Range[50];
P2=kfd1/.{x1->X[[#,1]],x2->X[[#,2]],x3->X[[#,3]]}&/@Range[51,100,1];

kfd1=kfd1-(Mean[P1]StandardDeviation[P2]+
Mean[P2]StandardDeviation[P1])/(StandardDeviation[P2]+
StandardDeviation[P1]);

kfdplot=ContourPlot3D[kfd1==0,{x1,2,8},{x2,0,6},{x3,0,6},
Mesh->None,ContourStyle->Opacity[0.3]];

Show[kfdplot,
ListPointPlot3D[{X[[1;;50]],X[[51;;100]]},
PlotStyle->PointSize[0.02],
PlotRange->{{2,8},{0,6},{0,6}}]
]

It's So Interesting!

Matehmatica 教學 Fisher Canonical Discriminant Analysis

2011-12-07T16:02:00.002+08:00

Fisher Canonical Discriminant Analysis

myline[x0_,y0_]:=
Line[{{x0,y0},{x,(-c-a*x)/b}/.
x->(b^2x0-a(c+b*y0))/(a^2+b^2)}]

data={RandomReal[
MultinormalDistribution[{1,1},{{2,1/2},{1/2,1}}],100],
RandomReal[MultinormalDistribution[{8,8},{{2,1/2},{1/2,1}}],100],
RandomReal[MultinormalDistribution[{3,8},{{2,1/2},{1/2,1}}],100]
};

Manipulate[
Block[{pts={point1,point2},a,b,c,m},
m=(pts[[2,2]]-pts[[1,2]])/(pts[[2,1]]-pts[[1,1]]);
b=-1;a=m;c=-m*pts[[1,1]]+pts[[1,2]];
Grid[{{a*x+b*y+c},
{ContourPlot[a*x+b*y+c==0,{x,-10,20},{y,-10,20},
ImageSize->350,
Epilog->{Red,Point@data[[1]],
Evaluate[myline@@@data[[1]]],
Green,Point@data[[2]],Evaluate[myline@@@data[[2]]],
Brown,Point@data[[3]],Evaluate[myline@@@data[[3]]]}
]}
}]
],{{point1,{-5,5}},Locator},{{point2,{15,15}},Locator}]

Mathematica 教學修改CSV格式並轉檔

2011-12-02T08:06:00.003+08:00

(*設定資料目錄*)
SetDirectory["~/Downloads"];
(*抓出所有的資料檔名*)
file = FileNames["*.txt"];
(*轉檔程式*)
recode[filename_String] := Block[{fn = filename, str},
(*讀入資料行，將文字轉成表達式，最後取出數字的元素*)
str = Cases[ToExpression@Import[fn, "Words"], _?NumberQ];
(*將資料依序輸出到xxx_new.txt檔*)
Do[str[[i]] >>> fn <> "_new.txt", {i, Length@str}]
]
(*將整個目錄的檔案轉檔*)
recode/@file

MS Word：論文中如何在『章』以中文編號，『節』以數字編號。

2011-11-17T15:27:00.003+08:00

今天上課時學生問了一個問題，如何在『章』以中文編號，『節』以數字編號。也就是說階層 1 用中文編號，階層 2, 階層 3, 階層 4 以後都以數字編號。

解決的方法很簡單，只要在『法律樣式編號』的地方打勾即可。見下圖：

Mathematica 教學 Path Analysis 路徑分析

2011-10-30T22:15:00.012+08:00

『懶』，是一種美德。
這個程式指定內生變數集外生變數後即可產生路徑圖以及路徑分析各項效果分解表，不需要在SPSS中作多次回歸或是在Amos中畫圖。

data=Import["http://www.pws.stu.edu.tw/cydye/Employee.csv"];

emp=data[[All,{2,4,5,6,7,8,9,10}]][[2;;-1]];

TableForm[emp[[1;;10]],
TableHeadings->{Range[Length@emp],
data[[1,{2,4,5,6,7,8,9,10}]]}]
path={{"educ","prevexp"},{"educ","salbegin"},{"educ",
"salary"},{"prevexp","salbegin"},{"prevexp",
"salary"},{"salbegin","salary"},{"jobtime",
"salary"},{"gender","prevexp"},{"gender",
"salbegin"},{"gender","salary"},{"minority",
"salary"},{"minority","salbegin"}};

exvars={"educ","gender","minority","jobtime"};
envars={"prevexp","salbegin","salary"};
vars=ToExpression@data[[1,{2,4,5,6,7,8,9,10}]];
MapThread[
Set,{vars,emp[[All,#]]&/@Range[Dimensions[emp][[-1]]]}];
pathtonode[targetnode_]:=Block[{node=targetnode,m,k},
NestWhile[DeleteCases[Flatten[Table[k=i;
m=Cases[path,{__,k[[1,1]]}];
Insert[k,#,1]&/@If[m=={},{{,}},m]
,{i,#}],1],{Null,Null},Infinity]&,{{{node,node}}},
Unequal,All][[All,1;;-2]]
]

pathtonode["salary"]

submodel[targetnode_]:=
Block[{node=targetnode,sourcenode,nodelist,varlist,pathcoeff},
sourcenode=
DeleteDuplicates[
Cases[pathtonode[ToString@node],{__,ToString@node},Infinity]];
nodelist=DeleteDuplicates[sourcenode][[All,1]];
nodelist=Insert[nodelist,node,-1];
varlist=
ToExpression@CharacterRange["a","z"][[1;;Length@nodelist-1]];
pathcoeff=
LinearModelFit[
Transpose[N@Standardize[#]&/@ToExpression@nodelist],
varlist,varlist]["ParameterTableEntries"][[All,{1,-1}]][[
2;;-1]];
MapThread[{Rule@@#1,
If[#2[[2]]<0.05, Style[NumberForm[#2[[1]],3, NumberFormat->(Row[{#1,
"\!\(\*SuperscriptBox[\(\),\(*\)]\)"}]&)],Bold,
FontSize->14],
Style[NumberForm[#2[[1]],3],Bold,FontSize->14]]}&,
{sourcenode,pathcoeff}]
]

pathgraph=
Flatten[DeleteCases[
If[pathtonode[#]!={{}},submodel[#]]&/@
DeleteDuplicates[Flatten[path]],Null],1];
Show[LayeredGraphPlot[pathgraph,Left,DirectedEdges->False],
LayeredGraphPlot[pathgraph,Left,VertexLabeling->True,
EdgeRenderingFunction->({Red,Arrow[#1,0.15]}&),
VertexRenderingFunction->({White,EdgeForm[Black],Black,
Text[Framed[Style[#2,FontSize->15],
Background->RGBColor[1,1,0.8]],#1]}&)],
ImageSize->350]


allpathrule=
Cases[pathgraph,{a_->b_,
Style[NumberForm[c_,3,___],Bold,FontSize->14]}:>{a,b}->
c,Infinity];

submodeleffect[sourcenode_,targetnode_]:=
Block[{tnode=targetnode,snode=sourcenode,tpathtonode,exeffect},
tpathtonode[ssourcenode_,ttargetnode_]:=
If[Length@Cases[envars,snode]==0,
Select[pathtonode[ttargetnode],#[[1,1]]==ssourcenode&],
DeleteDuplicates[
Flatten[{Cases[
pathtonode[ttargetnode],{___,{ssourcenode,__},___,{___,
ttargetnode},___},Infinity][[All,2;;-1]],
Cases[pathtonode[
ttargetnode],{___,{ssourcenode,ttargetnode},___},
Infinity][[All,-1]]},1]]];
exeffect=SortBy[tpathtonode[snode,tnode]/.allpathrule,Length@#&];
Which[Length@exeffect==0,Flatten@{0,0,0},
Length@exeffect==1,Flatten@{exeffect[[1]],0,exeffect[[1]]},
Length@exeffect>=2,Flatten@{exeffect[[1]],
Total[Times@@@exeffect[[2;;-1]]],
exeffect[[1]]+Total[Times@@@exeffect[[2;;-1]]]}]
]

Text@Grid[
Prepend[MapThread[
Insert[#1,#2,1]&,{Transpose[
Insert[Flatten[#]&/@
Partition[
Cases[Table[
submodeleffect[i,j],{j,envars},{i,
Flatten@{exvars,envars[[1;;-2]]}}],{x_,y_,
z_}:>{If[x>=0,""<>ToString@NumberForm[x,{4,3}],
NumberForm[x,{4,3}]],
If[y>=0,""<>ToString@NumberForm[y,{4,3}],
NumberForm[y,{4,3}]],
If[z>=0,""<>ToString@NumberForm[z,{4,3}],
NumberForm[z,{4,3}]]},Infinity],
Length@envars+Length@exvars-1],
Flatten[{"直接效果","間接效果","總效果"}&/@
Range[Length@envars+Length@exvars-1]],1]],
Flatten[{#,"",""}&/@Flatten[{exvars,envars[[1;;-2]]}]]}],
Flatten@{"","",envars}],
Background->{None,{Lighter[Yellow,0],{Lighter[Green,0.8],
Lighter[Blend[{Blue,Green}],.7],Lighter[Gray,0.5]}}},
Dividers->{{Darker[Gray,.6],{Lighter[Gray,.5]},
Darker[Gray,.6]},{Darker[Gray,.6],Darker[Gray,.6],{False},
Darker[Gray,.6]}},
Alignment->{Center},ItemSize->{{5,7,5,5,5}},
Frame->Darker[Gray,.6],ItemStyle->14,
Spacings->{Automatic,1}]

Finally, the second acceptance of 2011 comes.

2011-10-30T22:04:00.005+08:00

Application of particle swarm optimization for solving deteriorating inventory model with fluctuating demand and controllable deterioration rate

Yu-Ren Chen and Chung-Yuan Dye

Abstract

In most of the inventory models in the literature, the deterioration rate of goods is viewed as an exogenous variable, which is not subject to control. In the real market, the retailer can reduce the deterioration rate of product by making effective capital investment in storehouse equipments. In this study, we formulate a deteriorating inventory model with time-varying demand by allowing preservation technology cost as a deci- sion variable in conjunction with replacement policy. The objective is to find the optimal replenishment and preservation technology investment strategies while minimizing the total cost over the planning horizon. For any given feasible replenishment scheme, we first prove that the optimal preservation technology investment strategy not only exists but is unique. Then, a particle swarm optimization is coded and used to solve the non- linear programming problem by employing the properties derived from this paper. Some numerical examples are used to illustrate the features of the proposed model.

Key words: inventory, time-varying demand, deterioration, preservation technology investment, particle swarm optimization

Mathematica 教學更改Mathematica視窗設定

2011-10-26T07:27:00.007+08:00

SetOptions[InputNotebook[],
(* 視窗大小 *)
WindowSize -> {1200, 800},
(* 字型大小 *)
FontSize -> 18,
(* 輸出輸入是否要有外框 *)
CellFrame -> True,
(* 外框顏色 *)
CellFrameColor -> Green,
StyleDefinitions ->
Notebook[{Cell[StyleData[StyleDefinitions -> "Default.nb"]],
(* 輸入字型顏色 *)
Cell[StyleData["Input"], FontColor -> Blue],
(* 輸出字型顏色 *)
Cell[StyleData["Output"], FontColor -> Red]}]
]

若是要成為日後開啟的樣式，可將設定放在下面目錄的init.m裡面
ToFileName[$UserBaseDirectory, "Kernel"]

If[$Linked, (* do nothing for kernel-only sessions *)
RunScheduledTask[
FrontEndExecute[{
SetOptions[FrontEnd`InputNotebook[], WindowSize -> {1150,750}],
SetOptions[FrontEnd`InputNotebook[], WindowMargins ->
{{7,Automatic},{Automatic,6}}]
}], {0.4, 7}];
];

Finally, the first acceptance of 2011 comes.

2011-10-22T20:44:00.007+08:00

An optimal replenishment policy for deteriorating items with effective investment in preservation technology

Chung-Yuan Dye and Tsu-Pang Hsieh

Abstract

In this paper, considering the amount invested in preservation technology and the replenish- ment schedule as decision variables, we formulate an inventory model with a time-varying rate of deterioration and partial backlogging. The objective is to find the optimal replenishment and preservation technology investment strategies while maximizing the total profit per unit time. For any given preservation technology cost, we first prove that the optimal replenishment schedule not only exists but is unique. Next, under given replenishment schedule, we show that the total profit per unit time is a concave function of preservation technology cost. We then provide a simple algorithm to figure out the optimal preservation technology cost and replenishment schedule for the proposed model. We use numerical examples to illustrate the model.

Key words: inventory, deterioration, partial backlogging, preservation technology investment

Mathematica 教學 ContourPlot problem about color setting

2011-09-15T23:25:00.005+08:00

今天發現Mathematica的ContourPlot在處理一個比較複雜函數時若要求輸出ContourLabels時色彩的設定會出現問題

cp2 = ContourPlot[Evaluate@TP[t1, t2temp[t1, \[Xi]], \[Xi]], {t1, 0, 0.5}, {\[Xi], 0,
w}, Contours -> Range[17500, 18900, 200]]
cp3 = ContourPlot[Evaluate@TP[t1, t2temp[t1, \[Xi]], \[Xi]], {t1, 0, 0.5}, {\[Xi], 0,
w}, Contours -> Range[17500, 18900, 200], ContourLabels -> True]

比較以上兩個圖形會發現下面的圖形色彩的出現錯誤

(* 找出正常設定值 *)
zz = Cases[cp2[[1]], RGBColor[a_, b_, c_] :> {a, b, c}, Infinity];
z[i_] := zz[[i]]

(* 將正常設定值取代錯誤設定值 *)
(i = 1; Graphics[cp3[[1]] /. RGBColor[_, _, _] :> RGBColor @@ z[i++], cp3[[2]]])

Mathematica 教學 How to create a two-axis graph

2011-09-15T00:46:00.004+08:00

myplot[{f_,g_},range_List,color_List]:=
Block[{p1,p2,a1,a2,a3,b1,b2,myrange=range,mycolor=color},
p1=Plot[f,{x,myrange[[1]],myrange[[2]]},Axes->False,Frame->True];
p2=Plot[g,{x,myrange[[1]],myrange[[2]]},Axes->False,Frame->True];
a1=AbsoluteOptions[p1,FrameTicks][[1,2,2]];
a2=AbsoluteOptions[p2,FrameTicks][[1,2,2]];
b1=p1[[-1,4,-1,-1]];
b2=p2[[-1,4,-1,-1]];
a3=MapThread[
Insert[#1,#2,1]&,{a2[[All,2;;-1]],
b1[[1]]+((b1[[2]]-b1[[1]])/(b2[[2]]-b2[[1]])*(#-
b2[[1]]))&/@a2[[All,1]]}];

Plot[{f,b1[[1]]+((b1[[2]]-b1[[1]])/(b2[[2]]-b2[[1]])*
(g-b2[[1]]))},{x,myrange[[1]],myrange[[2]]},
PlotStyle->{{Thickness[0.01],mycolor[[1]]},
{Thickness[0.01],mycolor[[2]]}},
Axes->False,Frame->True,
FrameTicks->{{a1,a3},{Automatic,None}},
ImageSize->500]]

myplot[{0.25 Cos[x], 0.5 + Abs@Sin[x]}, {0, 2 Pi}, {Red, Blue}]

Mathematica 教學將矩陣依照指定規則重新組合

2011-09-05T23:34:00.003+08:00

Q: 依照矩陣元素位置轉換

Ans:


A = RandomInteger[100, {6, 6}]

transform[ma_List]:=Block[{myrule,tran},

(*轉換規則*)
myrule={{44,12,53,14,15,62},{21,22,23,43,25,26},
{31,13,33,34,35,36},{41,42,11,24,45,46},
{51,52,32,16,54,56},{61,55,63,64,65,66}};

(*取出對應位置元素*)
tran=ToExpression@
Flatten[Table[{StringTake[ToString[myrule[[i,j]]],1],
StringTake[ToString[myrule[[i,j]]],-1]},{i,6},{j,6}],1];

(*重新排列*)
Partition[ma[[#[[1]],#[[2]]]]&/@tran,6]]

(* test *)
transform[A] // MatrixForm

Mathematica 教學求解方程組後將解指定給各變數？

2011-09-01T09:07:00.006+08:00

這種方法只適用唯一解的情況！

myList={x,y};

ans=myList/.Solve[{x+y==a-b,2x-y==a+3b},myList][[1]];

MapThread[Set,{myList,ans}]

{x,y}