520表白嗨起来！高逼格5大统计软件勇表“脱单”心迹

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

论坛君：尽管笛卡尔和瑞典公主的故事已被证实只是杜撰，但因这个故事出名的心形函数被广为流传。今天又是一个虐单身狗的日子，面对各种毫无新意的表白方式，让我们来看看理工科式的表白~

---

一切都是从一个故事开始的：

1650年，斯德哥尔摩街头，一个宁静的午后，笛卡尔邂逅了18岁的瑞典公主克里斯汀。

机遇巧合，一段纯粹、美好的爱情悄然萌发。

然而，没过多久，他们的恋情传到了国王的耳朵里。国王大怒，将笛卡尔放逐，而公主被软禁在宫中。 身体孱弱的笛卡尔不久便染上重病，在生命进入倒计时的那段日子，他日夜思念的还是街头偶遇的那张温暖的笑脸。他每天坚持给她写信，盼望着她的回音。然而，这些信都被国王拦截下来。在笛卡尔给克里斯汀寄出第十三封信后，他永远地离开了这个世界。此时，被软禁在宫中的小公主依然徘徊在皇宫的走廊里，思念着远方的情人。

这最后一封信上没有写一句话，只有一个方程：r=a(1-sinθ) ——传说中著名的心形函数。

不过，事实告诉我们，除非你的目标妹子也是一只Geeker（至少会用Mathematica或者MATLAB等软件），否则像笛卡尔这样单给一个函数的结局大概就是你推公式别人推妹子了……表白什么的还是选择更浅显易懂的方法吧。

接下来看看用各种数据分析软件做出的心形模型，保证亮瞎你的表白对象！

Stata

在这个特殊的日子里，特意搜集了利用Stata绘制心形线的程序，代码如下：

1. clear
   
2. range t 0 2*_pi 1000
   
3. gen x=16*sin(t)^3
   
4. gen y=13*cos(t)-5*cos(2*t)-2*cos(3*t)-cos(4*t)
   
5. egen x_min=min(x)
   
6. egen x_max=max(x)
   
7. egen y_min=min(y)
   
8. egen y_max=max(y)
   
9. gen a=(x-x_min)/(x_max-x_min)
   
10. gen b=(y-y_min)/(y_max-y_min)
    
11. line b a
    
12. gr_edit yaxis1.draw_view.setstyle, style(no)
    
13. gr_edit xaxis1.draw_view.setstyle, style(no)
    
14. gr_edit plotregion1.AddTextBox added_text editor .7055394244311991 .2810707216715078
    
15. gr_edit plotregion1.added_text_new = 1
    
16. gr_edit plotregion1.added_text_rec = 1
    
17. gr_edit plotregion1.added_text[1].style.editstyle  angle(default) size(medsmall) color(red) horizontal(left) vertical(middle) margin(zero) linegap(zero) drawbox(no) boxmargin(zero) fillcolor(bluishgray) linestyle( width(thin) color(red) pattern(solid)) box_alignment(east) editcopy
    
18. gr_edit plotregion1.added_text[1].style.editstyle size(large) editcopy
    
19. gr_edit plotregion1.added_text[1].text = {}
    
20. gr_edit plotregion1.added_text[1].text.Arrpush  "           I LOVE YOU"
    
21. graph export "C:\Desktop\520.png", as(png) wid(800)hei(600) replace

复制代码

以下为输出结果：

以上程序基于stata12运行，如果有stata14的话，还可以输入中文！祝大家好运！

R

接下来是R绘制心形线，代码如下：

1. 1）加加载程序包
   
2. library(grid)
   
3. 
   
4. 2）绘制心形函数
   
5. heart <- function(lcolor){
   
6. t=seq(0, 2*pi, by=0.1)
   
7. x=16*sin(t)^3
   
8. y=13*cos(t)-5*cos(2*t)-2*cos(3*t)-cos(4*t)
   
9. a=(x-min(x))/(max(x)-min(x))
   
10. b=(y-min(y))/(max(y)-min(y))
    
11. grid.lines(a,b,gp=gpar(col=lcolor,lty = "solid",lwd = 6))
    
12. }
    
13. heart("hotpink")
    
14. grid.newpage()
    
15. 
    
16. 3）绘制玫瑰函数
    
17. rose=function(){
    
18. grid.circle(x=0.5, y=0.5, r=0.5,gp=gpar(fill="red",lwd = 3))
    
19. vp <- viewport(.5, .5, w=.9, h=.9)
    
20. pushViewport(vp)
    
21. grid.polygon(x=c(0.08, .5, 0.94),y=c(.22, 1.03, .22),gp=gpar(lwd =
    
22. vp2 <- viewport(.5, .5, w=.4, h=.4)
    
23. pushViewport(vp2)
    
24. grid.circle(x=0.5, y=0.5, r=0.5,gp=gpar(fill="red",lwd = 3))
    
25. vp3 <- viewport(.5, .5, w=.9, h=.9,angle=180)
    
26. pushViewport(vp3)
    
27. grid.polygon(x=c(0.08, .5, 0.94),y=c(.22, 1.03, .22),gp=gpar(lwd = 3))}
    
28. rose()
    
29. 
    
30. 4）调整图形
    
31. grid.newpage()
    
32. pushViewport(viewport(x=0.1, y=0.1,w=0.2, h=0.2))
    
33. grid.newpage()
    
34. for (j in 1:30) {
    
35. vp <- viewport(.5, .5, w=.9, h=.9)
    
36. pushViewport(vp)
    
37. heart("hotpink")
    
38. }
    
39. 
    
40. 5）调整图形及添加文字
    
41. grid.newpage()
    
42. vp1 <- viewport(.4, .5, w=.5, h=.5,angle=15)
    
43. pushViewport(vp1)
    
44. heart("red")
    
45. vp2 <- viewport(0.9, .27, w=.7, h=.7,angle=-30)
    
46. pushViewport(vp2)
    
47. heart("hotpink")
    
48. grid.text("执手偕老",
    
49.          x=0.2,y =1.2, just = c("center", "bottom"),
    
50.          gp = gpar(fontsize=30), vp = vp)
    
51. vp3 <- viewport(-0.65, 1.2, w=.3, h=.3,angle=-30)
    
52. pushViewport(vp3)
    
53. rose()

复制代码

结果如下：

注：操作平台为Rstudio，由于手头电脑没有装R，所以上述程序没有验证，就交给大家自己实践了！祝好运！

以上程序来自：
http://mp.weixin.qq.com/s?__biz=MzA4NTAyMjQ2Mg==&mid=206788501&idx=1&sn=53e07723fdb395a787e57ea53bffd332&scene=21#wechat_redirect

Matlab

利用matlab抱得美人归：

1. f=@(x,y,z) x.^2.*z.^3+9*y.^2.*z.^3/80-(x.^2+9*y.^2/4+z.^2-1).^3;%心形曲面函数
   
2. [x,y,z]=meshgrid(-1.5:0.1:1.5);%画图范围
   
3. v=f(x,y,z);
   
4. %画图
   
5. h=patch(isosurface(x,y,z,v,0));
   
6. isonormals(x,y,z,v,h)
   
7. set(h,'FaceColor','r','EdgeColor','none');
   
8. title('Programmed By Dylan Wang')
   
9. alpha(0.6)   
   
10. grid off;
    
11. axis([-1.5 1.5 -1.5 1.5 -1.5 1.5])
    
12. lighting Gouraud
    
13. h = camlight('left');
    
14. for i = 1:180;%水平旋转照相机
    
15.         camorbit(1,0)
    
16.         camlight(h,'left');
    
17.         drawnow;
    
18. end

复制代码

效果如斯：

SAS

1）准备心形图的图形数据

1. data heart;
   
2. length side $ 8;
   
3. pi=constant('PI');
   
4. do t=0 to 60 by 100/200;
   
5. x = -.01*(-t**2+40*t+1200)*sin(pi*t/180);
   
6. y = .01*(-t**2+40*t+1200)*cos(pi*t/180);
   
7. side='left';
   
8. output;
   
9. side='right';
   
10. x=-1*x;
    
11. output;
    
12. end;
    
13. run;

复制代码

2）准备图形上的文字

这个最重要，字一定要大，而且千万记得署上自己的名字。过年的时候收到几条没有署名的祝福短信，郁闷的我也不知道该怎么问人家名字。

1. data anno;
   
2. input textsize y1 function :$8. label $40.;
   
3. retain drawspace "GRAPHPERCENT"
   
4. function "text"
   
5. width 80
   
6. widthunit 'percent'
   
7. textcolor "black"
   
8. textweight "bold";
   
9. cards;
   
10. 20 50 text 5.20  我爱你！
    
11. 15 45 text 经管之家
    
12. ;
    
13. run;

复制代码

3）画图

1. ods listing;
   
2. ods graphics / width=600px height=600px noborder;
   
3. proc sgplot data=heart sganno=anno noautolegend;
   
4. title "Heart";
   
5. series x=x y=y /group=side lineattrs=(color=red);
   
6. xaxis display=none;
   
7. yaxis display=none;
   
8. run;quit;

复制代码

用SAS画曲线有多种方式：

1）用ScatterPlot画图是最自然的方式，不过需要10000个点心形图看上去才比较舒服。因为曲线的边有些锯齿，所以我不太喜欢。

2）用Series画图，这个需要点技巧。对数据不做特殊处理的话，画出来的是这样，不过也很美很特别。

Heart

3）我建议大家动手试试其他的语句，比如pbsplineplot，会有各种惊喜。

为了能出下面这个心形图，我特意加了side变量，然后通过lineattrs=(color=red)强行将两组数据用红色画线，而不是默认的一蓝一红。

Do you love SAS? Yes, I DO.（行动才是爱！）

Python

爱心的python表示：

1.第一种

1. """
   
2. '17*x^2 - 16*|x|*y + 17*y^2 = 225'
   
3. """
   
4. 
   
5. import numpy as np
   
6. import matplotlib.pyplot as plt
   
7. 
   
8. X = np.arange(-5.0, 5.0, 0.1)
   
9. Y = np.arange(-5.0, 5.0, 0.1)
   
10. 
    
11. x, y = np.meshgrid(X, Y)
    
12. f = 17 * x ** 2 - 16 * np.abs(x) * y + 17 * y ** 2 - 225
    
13. 
    
14. fig = plt.figure()
    
15. cs = plt.contour(x, y, f, 0, colors = 'r')
    
16. plt.show()

复制代码

2.第二种

1. """
   
2. '(x^2+y^2+y)^2 = x^2 + y^2'
   
3. """
   
4. import numpy as np
   
5. import matplotlib.pyplot as plt
   
6. 
   
7. X = np.arange(-2.0, 2.0, 0.05)
   
8. Y = np.arange(-2.0, 2.0, 0.05)
   
9. 
   
10. x, y = np.meshgrid(X, Y)
    
11. f = (x ** 2 + y ** 2 + y) ** 2 - x ** 2 - y ** 2
    
12. 
    
13. fig = plt.figure()
    
14. cs = plt.contour(x, y, f, 0, colors = 'r')
    
15. plt.show()

复制代码

3.第三种

1. """
   
2. '8*x^2 - 9*|x|*y + 8*y^2 = 17'
   
3. """
   
4. import numpy as np
   
5. import matplotlib.pyplot as plt
   
6. 
   
7. X = np.arange(-2.5, 2.5, 0.05)
   
8. Y = np.arange(-2.5, 2.5, 0.05)
   
9. 
   
10. x, y = np.meshgrid(X, Y)
    
11. f = 8 * x ** 2 - 9 * np.abs(x) * y + 8 * y ** 2 - 17
    
12. fig = plt.figure()
    
13. cs = plt.contour(x, y, f, 0, colors = 'r')

复制代码

4.第四种

1. """
   
2. '(x^2 + y^2 - 1)^3 - x^2*y^3 = 0'
   
3. """
   
4. import numpy as np
   
5. import matplotlib.pyplot as plt
   
6. import math
   
7. X = np.arange(-2.0, 2.0, 0.05)
   
8. Y = np.arange(-2.0, 2.0, 0.05)
   
9. 
   
10. x, y = np.meshgrid(X, Y)
    
11. 
    
12. f = (x ** 2 + y ** 2 - 1) ** 2 * (x ** 2 + y ** 2 - 1)- x ** 2 *  y ** 2 * y
    
13. fig = plt.figure()
    
14. cs = plt.contour(x, y, f, 0, colors = 'r')
    
15. plt.show()

复制代码

最后，感谢那些做出巨大贡献的程序猿，没有他们的代码，表白怎能如此风骚！！！


（关注微信+10个论坛币）
链接地址：https://bbs.pinggu.org/thread-2192087-1-1.html

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏7 回帖

关键词：统计软件 mathematica horizontal Mathematic Alignment 软件 统计

单身狗 路过！！！

来一个MATHEMATICA的呢，高手

ContourPlot3D[(x^2 + (9 y^2)/4 + z^2 - 1)^3 -
   x^2 z^3 - (9 y^2 z^3)/80 == 0, {x, -1.5, 1.5}, {y, -1.5,
  1.5}, {z, -1.5, 1.5}, PlotPoints -> 30, Axes -> False,
Lighting -> False, ContourStyle -> {RGBColor[1, 0, 0]}]

牛叉············

太有趣了

牛逼哄哄……

这个有点屌炸天