Stats

Popular Posts

Followers

這是留言板-既然來了,就留個言再走吧!

養花種魚數月亮賞星星 於 Monday, July 22, 2019 2:43 PM 發表


Please tell me what you thought about this blog. Comments, suggestions and views are welcomed.
繼續閱讀全文 這是留言板-既然來了,就留個言再走吧!

Mathematica 教學: Prove that 3/8< (1-cos(x)/cos(x/2))/x^2< 4/Pi^2 for all 0< x < Pi/2

Chung-Yuan Dye 於 Tuesday, September 19, 2017 8:52 PM 發表



KKTMinimize[obj_,eqns_,ineqns_,variables_]:=
Block[{myrule,lambda,u,\[Lambda],Lagrange,eqnu,ineqnlam,
eqnus,ineqnlams,kkteqns,kktvars,kktans},
myrule={z_[i_]:>ToExpression[ToString[z]<>ToString[i]]};
eqnu=u[#]&/@Range[Length@eqns];
ineqnlam=\[Lambda][#]&/@Range[Length@ineqns];
eqnus=If[Length@eqns>=1,eqns.eqnu,0];
ineqnlams=If[Length@ineqns>=1,ineqns.ineqnlam,0];
kktvars=Flatten@{variables,eqnu,ineqnlam}/.myrule;
Lagrange=obj-eqnus-ineqnlams;
kkteqns=Flatten@{Thread[D[Lagrange,{variables}]==0],
Thread[eqns==0],Thread[ineqns<=0],
Thread[ineqns*ineqnlam==0],Thread[ineqnlam<=0]}/.
myrule;
If[MemberQ[PolynomialQ[#,kktvars]&/@kkteqns[[All,1]],False],
Print["KKT限制式均需為多項式。"],
kktans=Reduce[kkteqns,kktvars,
Backsubstitution->True]/.{And->List,Or->List,
Equal->Rule};
If[Length@Dimensions@kktans==1,{obj/.kktans,
kktans},{obj/.#,#}&/@kktans]]];

Plot[(1-Cos[x]/Cos[x/2])/x^2,{x,0.-10,Pi/2},
PlotRange->{{0,Pi/2},{0.3,0.5}},Frame->True,
FrameTicks->{{{0.3,0.5,{3/8,"3/8"},{4/Pi^2,"4/Pi^2"}},
None},{{0,Pi/2},None}},
GridLines->{None,{3/8,4/Pi^2}},GridLinesStyle->Dashed]

D[(1-Cos[x]/Cos[x/2])/x^2,x]//Simplify

D[(1-(Cos[x/2]^2-Sin[x/2]^2)/Cos[x/2])/x^2,x]//Simplify

D[4-4Sec[x/2]+3xTan[x/2]-4Tan[x/2]^2+xTan[x/2]^3,x]//TrigExpand//Simplify

(D[(6x-2Sin[x/2]-4Sin[x]-2Sin[(3x)/2]+Sin[2x]),x]//TrigExpand)
/.{Sin[z_]:>a,Cos[z_]:>b}

KKTMinimize[6+4a^2+2a^4-b+9a^2b-4b^2-12a^2b^2-3b^3+2b^4,{},
{-a,-b,a-1/Sqrt[2],1/Sqrt[2]-b,b-1},{a,b}]

Limit[(1-(Cos[x/2]^2-Sin[x/2]^2)/Cos[x/2])/x^2,x->0]

Limit[(1-(Cos[x/2]^2-Sin[x/2]^2)/Cos[x/2])/x^2,x->Pi/2]

Let f(x)=(1-cos(x)/cos(x/2))/x^2. Taking the first derivative of f(x) with respect to x yields d/dx f(x)=cos(x/2)*g(x)/(2x^3),
where g(x)=x*tan^3(x/2)-4*tan^2(x/2)+3*x*tan(x/2)-4*sec(x/2)+4. Because g'(x)=1/4*sec(x/2)^2*h(x) and h'(x)=k(x)=2*sin^4(x/2)+4*sin^2(x/2)+2*cos^4(x/2)-3*cos^3(x/2)-4*cos^2(x/2)-cos(x/2)-12*sin^2(x/2) *cos^2(x/2)+9*sin^2(x/2) cos(x/2)+6. Let a=sin(x/2) and b=cos(x/2). By Karush-Kuhn-Tucker conditions, we know that k(x) reaches its minimum at a=0 and b=1. This, together with the facts that h(0)=0, g(0)=0, and f(0)=0, indicates that f(x) is increasing in x, which also implies that 3/8<(1-Cos[x]/Cos[x/2])/x^2<4/π^2. This completes the proof.
繼續閱讀全文 Mathematica 教學: Prove that 3/8< (1-cos(x)/cos(x/2))/x^2< 4/Pi^2 for all 0< x < Pi/2

Mathematica 教學: Prove that 1 - 3 s^2 + (1 + s^2) Cosh[s] >0.

Chung-Yuan Dye 於 Monday, September 18, 2017 11:55 PM 發表



Plot[1-3s^2+(1+s^2)Cosh[s],{s,0,Pi},Frame->True]
SeriesCoefficient[1-3s^2+(1+s^2)Cosh[s],{s,0,n}]
Normal[Series[1-3s^2+(1+s^2)Cosh[s],{s,0,5}]]
Minimize[2-(3s^2)/2+(13s^4)/24,s]
Solve[D[2-(3s^2)/2+(13s^4)/24,s]==0,s]

Since the coefficient in the n-th order of the taylor series of 1-3s^2+(1 + s^2)Cosh[s] is positive for every even term and is
zero for every odd term for all n>=3, it is obvious to see that 1-3s^2+(1+s^2)Cosh[s]>2-(3s^2)/2+(13s^4)/24>25/26, which implies that 1-3s^2+(1+s^2)Cosh[s]>0. This completes the proof.
繼續閱讀全文 Mathematica 教學: Prove that 1 - 3 s^2 + (1 + s^2) Cosh[s] >0.

Mathematica 教學:抓取證交所股價資料

Chung-Yuan Dye 於 Sunday, July 23, 2017 11:36 PM 發表


mydate[date_]:=Block[{temp},
temp=Import[
"http://www.twse.com.tw/exchangeReport/MI_INDEX.html?date="<>ToString@date<>"&type=ALLBUT0999",
"JSON"];
{("data1"/.temp)[[All,{1,2,4,5}]],("data5"/.temp)[[All,
{1,2,3,4,5,6,7,8,9,11}]]}]

mydate[20170719][[1]] // TableForm
mydate[20170719][[2]] // TableForm

test=Quiet@{
{ToExpression@StringJoin@Characters[ToString@#][[1;;4]],
ToExpression@StringJoin@Characters[ToString@#][[5;;6]],
ToExpression@StringJoin@Characters[ToString@#][[7;;8]]},
mydate[#][[1,4]]}&/@Map[ToExpression@
StringJoin[ToString/@{#[[1]],
If[#[[2]]<10>ToString[#[[2]]],ToString[#[[2]]]],
If[#[[3]]<10>ToString[#[[3]]],
ToString[#[[3]]]]}]&,
DayRange["2017/05/1","2017/07/21"][[All,1]]]//Parallelize;

DateListPlot[{#[[1]],
ToExpression@StringReplace[#[[2,2]],","->""]}&/@
Select[test,Length@#[[2]]==4&],Filling->Bottom]


Block[{tempdata},
tempdata={#[[1]],ToExpression@StringReplace[#[[2,2]],","->""]}&/@
Select[test,Length@#[[2]]==4&];
DateListPlot[tempdata,Filling->Bottom,
Epilog->{Red,PointSize[0.025],Point[{tempdata[[#[[1]]]][[1]],#[[2]]}]&/@
FindPeaks[tempdata[[All,2]]]}]]
 

繼續閱讀全文 Mathematica 教學:抓取證交所股價資料

ListPlot with Colormap

Chung-Yuan Dye 於 Friday, June 30, 2017 8:26 PM 發表

mylistplot[data__,xrange_,yrange_,pointsize_,colormap_]:=
Block[{temp},
temp=ListPlot[data,
PlotRange->{xrange,yrange},
PlotStyle->PointSize[pointsize],
Frame->True,
ColorFunction->colormap,
Joined->True][[1,2,3,2]];

ListPlot[
Style[#[[1]],RGBColor[#[[2]]]]&/@Transpose[{data,temp}],
PlotTheme->"Detailed",
PlotStyle->PointSize[pointsize],
PlotRange->{xrange,yrange},
Frame->True,
PlotLegends->Placed[BarLegend[{colormap,yrange}],Right]]
]

mylistplot[iris,{2,4},{4,8},0.015,"TemperatureMap"]

繼續閱讀全文 ListPlot with Colormap

Use R within Mathematica

Chung-Yuan Dye 於 Thursday, June 22, 2017 10:22 PM 發表

Needs["RLink`"]

(* For Mac and MMA 10*)

SetEnvironment[ "DYLD_LIBRARY_PATH" ->"/Library/Frameworks/R.framework/Resources/lib"]
InstallR["RHomeLocation" -> "/Library/Frameworks/R.framework/Resources", "RVersion" -> 3];
REvaluate["R.version.string"]

(* Load packages*)

REvaluate["library(psych);"]
REvaluate["library(nortest);"]
REvaluate["library(GenABEL);"]

(* 因素分析,主成分法萃取,變異及大化法轉軸 *)

myfac=Map[REvaluate["{xx<-principal(iris[,1:4],nfactors="<>ToString[#]<>",
rotate=\"varimax\");xx$loadings}"]&,{2,3,4}];
myfac

(* 輸出因素負荷矩陣 *)

Reverse[SortBy[myfac[[1,1]],{#[[1]],#[[2]]}&]]//TableForm
TableForm/@myfac[[All,1]]


(* 因素分析資料*) 

fadata = {
{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 4, 5, 6}, 
{1, 2, 1, 1, 1, 1, 2, 1, 2, 1, 3, 4, 3, 3, 3, 4, 6, 5}, 
{3, 3, 3, 3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 5, 4, 6}, 
{3, 3, 4, 3, 3, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 5, 6, 4}, 
{1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1, 6, 4, 5}, 
{1, 1, 1, 2, 1, 3, 3, 3, 4, 3, 1, 1, 1, 2, 1, 6, 5, 4}};

(* 定義資料*)

RSet["factor", Transpose@fadata];

(* 在R中計算相關係數矩陣,傳回MMA小數後四捨五入 *)

Round[REvaluate["cor(factor)"], 0.0001] // MatrixForm

(* 因素分析 *)

REvaluate["factanal(factor, factors = 3)$loadings"]

(* 輸出因素負荷矩陣 *)

TableForm[#[[1]], TableHeadings -> {None, #[[-1, -1, -1, -1]]}] & /@ 
{REvaluate["factanal(factor, factors = 3)$loadings"]} // Column
 
 
(* 迴歸分析 *)

REvaluate["{
data(iris)
reg <- lm( Sepal.Length ~ Species, data=iris )
summary.text <- capture.output(print( summary(reg)) )
}"]

StringJoin@Riffle[#, "\n"] &@REvaluate["{
data(iris)
reg <- lm( Sepal.Length ~ Species, data=iris )
summary.text <- capture.output(print( summary(reg)) )
}"]
 
(* 產生常態隨機變數 *)

mydata=REvaluate["rnorm(100,mean=0,sd=1)"]&/@Range[100]; 

(* MMA次數分配 *)

\[ScriptCapitalD] = HistogramDistribution[mydata[[1]]]  

(* MMA PDF & CDF *)

GraphicsRow[DiscretePlot[#[\[ScriptCapitalD], x], {x, -4, 4, .01}, 
PlotLabel -> #] & /@ {PDF, CDF}]

(* shapiro test, ks.test *)

shapiro[data_] := Block[{temp},
  RSet["temp", data];
  REvaluate[
 "{xx=shapiro.test(temp); yy=ks.test(temp,\"pnorm\",0,1); zz=lillie.test(temp); 
 c(xx$statistic,xx$p.value,yy$statistic,yy$p.value,zz$statistic,zz$p.value)}"][[1]] ]
 
RSet["rn1", RandomReal[{0, 10}, 1000]];
rn2 = RandomReal[NormalDistribution[0, 1], 1000];
REvaluate["rntransform(rn1)"] // Short
Histogram[#, Frame -> True] & /@ {rn2, REvaluate["rntransform(rn1)"]}
ShapiroWilkTest[#, {"TestStatistic", "PValue"}] & /@ {rn2, REvaluate["rntransform(rn1)"]} 


(* 利用RCurl擷取網路資料 *)

gosspost[n_]:=REvaluate["
{library(RCurl);
curl<-getCurlHandle();
curlSetOpt(cookie=\"over18=1\",followlocation=TRUE,curl=curl);
url<-\"https://www.ptt.cc/bbs/Gossiping/index" <> ToString[n] <> 
".html\";
html<-getURL(url,curl=curl)
}"][[1]];

Gossiping = 
  Cases[ImportString[
StringReplace[gosspost[#], {"\n" -> "", "\t" -> ""}], 
"XMLObject"],
  XMLElement[
 "div", {"class" -> "r-ent"}, {XMLElement[
"div", {"class" -> "nrec"}, {XMLElement[
  "span", {"class" -> "hl f2"}, {_}]}], 
  XMLElement["div", {"class" -> "mark"}, {}], 
  XMLElement[
"div", {"class" -> "title"}, {XMLElement[
  "a", {"shape" -> "rect", 
"href" -> postlink_}, {posttile_}]}], 
  XMLElement[
"div", {"class" -> "meta"}, {XMLElement[
  "div", {"class" -> "date"}, {_}], 
 XMLElement[
  "div", {"class" -> "author"}, {author_}]}]}] :> {postlink, 
 posttile, author},
  Infinity] & /@ Range[20861, 20862, 1];
  
Flatten[Gossiping, 1] // TableForm  


繼續閱讀全文 Use R within Mathematica

Thesis No Worry

Chung-Yuan Dye 於 Thursday, April 6, 2017 2:11 PM 發表


thesisnoworryhlm[構面1量表,構面2量表,構面3量表,{構面1各分量表題數},{構面2各分量表題數},{構面3各分量表題數},
{構面1各分量表變數名稱},{構面2各分量表變數名稱},{構面3各分量表變數名稱}]

thesisnoworryhlm[ddata[[All,1;;15]],ddata[[All,16;;24]],
ddata[[All,25;;30]],{5,6,4},{4,5},{3,3},
{X1,X2,X3},{M1,M2},{Y1,Y2}]


thesisnoworryanova[構面1量表,構面2量表,構面3量表,{構面1各分量表題數},{構面2各分量表題數},{構面3各分量表題數},
{構面1各分量表變數名稱},{構面2各分量表變數名稱},{構面3各分量表變數名稱},人口統計變量資料,{人口統計變量變數名稱}


thesisnoworryanova[ddata[[All,1;;15]],ddata[[All,16;;24]],
ddata[[All,25;;30]],{5,6,4},{4,5},{3,3},
{X1,X2,X3},{M1,M2},{Y1,Y2},
ddata[[All,{31,32,33}]],{"Sex","Age","Edu"}]



Thesis No Worry 程式資料檔

繼續閱讀全文 Thesis No Worry

How to get the coefficient of polynomial function

Chung-Yuan Dye 於 Monday, October 3, 2016 6:52 PM 發表


1+2Sin[2x]+3Sin[3x]+4Sin[4x]/.{Plus->List,Sin[y_]->1}
繼續閱讀全文 How to get the coefficient of polynomial function

原生植物園

Chung-Yuan Dye 於 Sunday, September 18, 2016 9:11 PM 發表






















繼續閱讀全文 原生植物園

落霞與孤鶩齊飛,秋水共長天一色。

Chung-Yuan Dye 於 Saturday, September 17, 2016 12:18 AM 發表
落霞與孤鶩齊飛,秋水共長天一色。可惜這裡不是鄱陽湖的滕王閣,而是蓮池潭的五里亭!


















繼續閱讀全文 落霞與孤鶩齊飛,秋水共長天一色。

樹德隨拍

Chung-Yuan Dye 於 Friday, September 9, 2016 11:45 PM 發表



















繼續閱讀全文 樹德隨拍

雨後隨拍

Chung-Yuan Dye 於 Thursday, September 8, 2016 3:39 PM 發表

繼續閱讀全文 雨後隨拍

樹德科技大學寄情胡

Chung-Yuan Dye 於 Monday, August 29, 2016 4:47 PM 發表










繼續閱讀全文 樹德科技大學寄情胡