График распределения Гаусса в 3D

10

В теории вероятностей нормальное (или гауссово) распределение является очень распространенным непрерывным распределением вероятности. Нормальные распределения важны в статистике и часто используются в естественных и социальных науках для представления вещественных случайных величин, распределение которых неизвестно.

Соревнование

Ваша задача - построить плотность вероятности распределения Гаусса на трехмерной плоскости . Эта функция определяется как:

Куда:




A = 1, σ x = σ y = σ

правила

  • Ваша программа должна принимать один вход σ , стандартное отклонение.
  • Ваша программа должна распечатать трехмерный график распределения Гаусса в максимально высоком качестве, насколько позволяет ваш язык / система.
  • Ваша программа не может использовать встроенное прямое распределение Гаусса или плотность вероятности.
  • Ваша программа не должна прекращаться.
  • Ваш сюжет может быть черно-белым или цветным.
  • Ваш график должен иметь линии сетки внизу. Линии сетки по бокам (как показано в примерах) не нужны.
  • Ваш график не должен иметь номера строк рядом с линиями сетки.

счет

Как обычно в , выигрывает представление с наименьшим количеством байтов! Я, возможно, никогда не «приму» ответ с помощью кнопки, если только он невероятно мал и интуитивен.

Пример вывода

Ваш вывод может выглядеть примерно так:

5

Или это может выглядеть так:

6

Более действительные выводы . Неверные выводы .

MD XF
источник
Я был смущен, что вы только что показали функцию для оси X. Нужно ли брать отдельные входы / выходы для X и Y сигмы и мю?
Скотт Милнер
Таким образом, мы должны предполагать, что μ равно 0? А какой масштаб вам нужен для х и у? Если x- и y-диапазоны выбраны очень малыми относительно σ, то график по существу будет выглядеть как постоянная функция.
Грег Мартин
(Я думаю, что для двумерного распределения будет понятнее, если вы используете | x-μ | ^ 2 в определении, а не (x-μ) ^ 2.)
Грег Мартин
@GregMartin Отредактировано.
MD XF
2
Все еще не ясно ... что такое x_o, y_o и θ?
Грег Мартин

Ответы:

7

Гнуплот 4, 64 62 61 60 47 байт

(Связано с Mathematica ! WooHoo!)

se t pn;se is 80;sp exp(-(x**2+y**2)/(2*$0**2))

Сохраните приведенный выше код в файл с именем A.gpи вызовите его со следующим:

gnuplot -e 'call "A.gp" $1'>GnuPlot3D.png

где $1должен быть заменен значением σ. Это сохранит .pngфайл с именем, GnuPlot3D.pngсодержащим желаемый вывод, в текущий рабочий каталог.

Обратите внимание, что это работает только с дистрибутивами Gnuplot 4, поскольку в Gnuplot 5 $nссылки на аргументы устарели и заменены, к сожалению, более многословными ARGn.

Пример вывода с σ = 3:

Пример вывода

Этот вывод в порядке в соответствии с OP .


Гнуплот 4, Альтернативное решение, 60 байтов

Вот альтернативное решение, которое намного длиннее предыдущего, но на мой взгляд вывод выглядит намного лучше.

se t pn;se is 80;se xyp 0;sp exp(-(x**2+y**2)/(2*$0**2))w pm

Это все еще требует Gnuplot 4 по той же причине, что и предыдущее решение.

Пример вывода с σ = 3:

Пример вывода № 2

Р. Кап
источник
I am not sure if it molds to the specifications requiredкакие технические характеристики, по вашему мнению, не соответствуют?
MD XF
@MDXF Во-первых, я не уверен, что прозрачность графика в порядке. Честно говоря, мне это не очень нравится, поэтому я не был уверен, что здесь все будет в порядке. Во-вторых, по умолчанию график начинается на одну единицу сверху вниз, и я не уверен, что все в порядке. В-третьих, поскольку график начинается на одну единицу выше, я не уверен, что диспропорциональность графика по сравнению с графиками, приведенными в оригинальном посте, в порядке. Однако, если с тобой все в порядке, я с радостью сделаю это главным ответом.
Р. Кап
@MDXF На самом деле я собирался опубликовать его как исходный ответ, но по этим причинам я предпочел не делать этого, а вместо этого опубликовал текущий ответ.
Р. Кап
@MDXF На самом деле, я могу сделать еще короче , если это хорошо. Я понимаю, если это не будет, но это не больно спрашивать. Это способ по умолчанию Gnuplotпостроить график плотности вероятности распределения Гаусса с сигма 2без каких-либо модификаций среды.
Р. Кап
@MDXF Думаю, я мог бы спросить, прежде чем публиковать свой первоначальный ответ, но в то время я очень хотел опубликовать ответ.
Р. Кап
14

C ++, 3477 3344 байта

Количество байтов не включает ненужные переводы строки.
MD XF забил 133 байта.

Нет никакого способа, которым C ++ мог бы конкурировать за это, но я думал, что было бы забавно написать программное средство рендеринга для задачи. Я разорвал и сыграл несколько частей GLM для математики 3D и использовал линейный алгоритм Сяолиня Ву для растеризации. Программа выводит результат в файл PGM с именем g.

Вывод

#include<array>
#include<cmath>
#include<vector>
#include<string>
#include<fstream>
#include<algorithm>
#include<functional>
#define L for
#define A auto
#define E swap
#define F float
#define U using
U namespace std;
#define K vector
#define N <<"\n"
#define Z size_t
#define R return
#define B uint8_t
#define I uint32_t
#define P operator
#define W(V)<<V<<' '
#define Y template<Z C>
#define G(O)Y vc<C>P O(vc<C>v,F s){vc<C>o;L(Z i=0;i<C;++i){o\
[i]=v[i]O s;}R o;}Y vc<C>P O(vc<C>l, vc<C>r){vc<C>o;L(Z i=0;i<C;++i){o[i]=l[i]O r[i];}R o;}
Y U vc=array<F,C>;U v2=vc<2>;U v3=vc<3>;U v4=vc<4>;U m4=array<v4,4>;G(+)G(-)G(*)G(/)Y F d(
vc<C>a,vc<C>b){F o=0;L(Z i=0;i<C;++i){o+=a[i]*b[i];}R o;}Y vc<C>n(vc<C>v){R v/sqrt(d(v,v));
}v3 cr(v3 a,v3 b){R v3{a[1]*b[2]-b[1]*a[2],a[2]*b[0]-b[2]*a[0],a[0]*b[1]-b[0]*a[1]};}m4 P*(
m4 l,m4 r){R{l[0]*r[0][0]+l[1]*r[0][1]+l[2]*r[0][2]+l[3]*r[0][3],l[0]*r[1][0]+l[1]*r[1][1]+
l[2]*r[1][2]+l[3]*r[1][3],l[0]*r[2][0]+l[1]*r[2][1]+l[2]*r[2][2]+l[3]*r[2][3],l[0]*r[3][0]+
l[1]*r[3][1]+l[2]*r[3][2]+l[3]*r[3][3]};}v4 P*(m4 m,v4 v){R v4{m[0][0]*v[0]+m[1][0]*v[1]+m[
2][0]*v[2]+m[3][0]*v[3],m[0][1]*v[0]+m[1][1]*v[1]+m[2][1]*v[2]+m[3][1]*v[3],m[0][2]*v[0]+m[
1][2]*v[1]+m[2][2]*v[2]+m[3][2]*v[3],m[0][3]*v[0]+m[1][3]*v[1]+m[2][3]*v[2]+m[3][3]*v[3]};}
m4 at(v3 a,v3 b,v3 c){A f=n(b-a);A s=n(cr(f,c));A u=cr(s,f);A o=m4{1,0,0,0,0,1,0,0,0,0,1,0,
0,0,0,1};o[0][0]=s[0];o[1][0]=s[1];o[2][0]=s[2];o[0][1]=u[0];o[1][1]=u[1];o[2][1]=u[2];o[0]
[2]=-f[0];o[1][2]=-f[1];o[2][2]=-f[2];o[3][0]=-d(s,a);o[3][1]=-d(u,a);o[3][2]=d(f,a);R o;}
m4 pr(F f,F a,F b,F c){F t=tan(f*.5f);m4 o{};o[0][0]=1.f/(t*a);o[1][1]=1.f/t;o[2][3]=-1;o[2
][2]=c/(b-c);o[3][2]=-(c*b)/(c-b);R o;}F lr(F a,F b,F t){R fma(t,b,fma(-t,a,a));}F fp(F f){
R f<0?1-(f-floor(f)):f-floor(f);}F rf(F f){R 1-fp(f);}struct S{I w,h; K<F> f;S(I w,I h):w{w
},h{h},f(w*h){}F&P[](pair<I,I>c){static F z;z=0;Z i=c.first*w+c.second;R i<f.size()?f[i]:z;
}F*b(){R f.data();}Y vc<C>n(vc<C>v){v[0]=lr((F)w*.5f,(F)w,v[0]);v[1]=lr((F)h*.5f,(F)h,-v[1]
);R v;}};I xe(S&f,v2 v,bool s,F g,F c,F*q=0){I p=(I)round(v[0]);A ye=v[1]+g*(p-v[0]);A xd=
rf(v[0]+.5f);A x=p;A y=(I)ye;(s?f[{y,x}]:f[{x,y}])+=(rf(ye)*xd)*c;(s?f[{y+1,x}]:f[{x,y+1}])
+=(fp(ye)*xd)*c;if(q){*q=ye+g;}R x;}K<v4> g(F i,I r,function<v4(F,F)>f){K<v4>g;F p=i*.5f;F
q=1.f/r;L(Z zi=0;zi<r;++zi){F z=lr(-p,p,zi*q);L(Z h=0;h<r;++h){F x=lr(-p,p,h*q);g.push_back
(f(x,z));}}R g;}B xw(S&f,v2 b,v2 e,F c){E(b[0],b[1]);E(e[0],e[1]);A s=abs(e[1]-b[1])>abs
(e[0]-b[0]);if(s){E(b[0],b[1]);E(e[0],e[1]);}if(b[0]>e[0]){E(b[0],e[0]);E(b[1],e[1]);}F yi=
0;A d=e-b;A g=d[0]?d[1]/d[0]:1;A xB=xe(f,b,s,g,c,&yi);A xE=xe(f,e,s,g,c);L(I x=xB+1;x<xE;++
x){(s?f[{(I)yi,x}]:f[{x,(I)yi}])+=rf(yi)*c;(s?f[{(I)yi+1,x}]:f[{x,(I)yi+1}])+=fp(yi)*c;yi+=
g;}}v4 tp(S&s,m4 m,v4 v){v=m*v;R s.n(v/v[3]);}main(){F l=6;Z c=64;A J=g(l,c,[](F x,F z){R
v4{x,exp(-(pow(x,2)+pow(z,2))/(2*pow(0.75f,2))),z,1};});I w=1024;I h=w;S s(w,h);m4 m=pr(
1.0472f,(F)w/(F)h,3.5f,11.4f)*at({4.8f,3,4.8f},{0,0,0},{0,1,0});L(Z j=0;j<c;++j){L(Z i=0;i<
c;++i){Z id=j*c+i;A p=tp(s,m,J[id]);A dp=[&](Z o){A e=tp(s,m,J[id+o]);F v=(p[2]+e[2])*0.5f;
xw(s,{p[0],p[1]},{e[0],e[1]},1.f-v);};if(i<c-1){dp(1);}if(j<c-1){dp(c);}}}K<B> b(w*h);L(Z i
=0;i<b.size();++i){b[i]=(B)round((1-min(max(s.b()[i],0.f),1.f))*255);}ofstream f("g");f 
W("P2")N;f W(w)W(h)N;f W(255)N;L(I y=0;y<h;++y){L(I x=0;x<w;++x)f W((I)b[y*w+x]);f N;}R 0;}
  • l длина одной стороны сетки в мировом пространстве.
  • c это число вершин вдоль каждого края сетки.
  • Функция, которая создает сетку, вызывается с функцией, которая принимает два входа, xи z(+ y идет вверх), мировые пространственные координаты вершины и возвращает мировую пространственную позицию вершины.
  • w это ширина пгм
  • h высота пгм
  • mматрица вида / проекции Аргументы, используемые для создания m...
    • поле зрения в радианах
    • соотношение сторон pgm
    • ближняя плоскость клипа
    • дальний отсек плоскости
    • положение камеры
    • цель камеры
    • вверх вектор

Рендерер может легко иметь больше возможностей, лучшую производительность и лучше играть в гольф, но мне было весело!

Патрик Перселл
источник
2
Вау, это невероятно!
MD XF
1
Вовсе нет ... дерзайте!
Патрик Перселл,
1
Вот и все, 133 байта!
MD XF
1
Это потрясающе! Если бы вы могли сказать мне, где вы узнали все это, это было бы здорово !
HatsuPointerKun
1
@HatsuPointerKun Рад, что вам понравится! Этот учебник ... opengl-tutorial.org/beginners-tutorials/tutorial-3-matrices ... является отличным местом для начала.
Патрик Перселл,
9

Mathematica, 47 байт

Plot3D[E^(-(x^2+y^2)/2/#^2),{x,-6,6},{y,-6,6}]&

принимает в качестве входа σ

вход

[2]

вывод
введите описание изображения здесь

-2 байта благодаря LLlAMnYP

J42161217
источник
1
Mathematica выигрывает? Там нет сюрпризов: P
MD XF
3
Сохранение 2 байтов с помощьюE^(-(x^2+y^2)/2/#^2)
LLlAMnYP
6

R 105 102 87 86 байт

s=scan();plot3D::persp3D(z=sapply(x<-seq(-6,6,.1),function(y)exp(-(y^2+x^2)/(2*s^2))))

Принимает сигму от STDIN. Создает вектор из -6в 6в шагах .1для обоих xи y, затем создает 121x121матрицу, беря внешнее произведение xи y. Это короче, чем вызов matrixи указание размеров. Матрица уже заполнена, но это нормально, потому что мы перезаписываем это.

В for-loop петли над значениями в x, что делает использование векторизованных операций R, создающая матрицу плотности по одной строке за один раз.

(s)applyснова более короткий метод для векторизованных операций. Как и герой, он сам обрабатывает создание матрицы, сохраняя немало байтов.

введите описание изображения здесь

128 125 110 109 байт, но гораздо более навороченный:

Этот сюжет создан plotlyпакетом. К сожалению, спецификация немного многословна, так что это стоит много байтов. Результат действительно очень интересный. Я очень рекомендую попробовать это для себя.

s=scan();plotly::plot_ly(z=sapply(x<-seq(-6,6,.1),function(y)exp(-(y^2+x^2)/(2*s^2))),x=x,y=x,type="surface")

бла

JAD
источник
Я указал в вопросе, что график не должен иметь номера строк, ваше второе представление в порядке.
MD XF
О, должно быть, я пропустил это. Я обменял свои решения вокруг. Я думаю, что plotlyсюжет достаточно причудлив, чтобы оправдать его включение.
JAD
Ну, и то, и другое намного, намного красивее моего : P
MD XF
Поскольку вы используете только sодин раз, вы могли бы сделать 2*scan()^2и удалить s=scan();в начале? Это сэкономит 3 байта.
KSmarts
6

Applesoft BASIC, 930 783 782 727 719 702 695 637 байт

-72 байта и работающая программа благодаря потолку, определяющему мою ошибку и укороченному алгоритму

0TEXT:HOME:INPUTN:HGR:HCOLOR=3:W=279:H=159:L=W-100:Z=L/10:B=H-100:C=H-60:K=0.5:M=1/(2*3.14159265*N*N):FORI=0TO10STEPK:X=10*I+1:Y=10*I+B:HPLOTX,Y:FORJ=0TOL STEP1:O=10*J/L:D=ABS(5-I):E=ABS(5-O):R=(D*D+E*E)/(2*N*N):G=EXP(-R)*M:A=INT((C*G)/M):X=10*I+Z*O+1:Y=10*I+B-A:HPLOTTOX,Y:IF(I=0)GOTO4
1IF(J=L)GOTO3
2V=INT(J/10):IF((J/10)<>V)GOTO5
3D=ABS(5-I+K):E=ABS(5-O):R=(D*D+E*E)/(2*N*N):U=EXP(-R)/(2*3.14159*N*N):S=INT((C*U)/M):P=10*(I-K)+Z*O+1:Q=10*(I-K)+B-S:HPLOT TOP,Q:HPLOTX,Y
4IF(J=0)GOTO7:IF(I<10)GOTO5:IF(J=L)GOTO6:V=INT(J/10):IF((J/10)=V)GOTO6
5HCOLOR=0
6HPLOTTOX,10*I+B:HCOLOR=3:HPLOTX,Y
7NEXTJ:NEXTI:HPLOTW+1,H:HPLOTTO101,H:HPLOTTO0+1,H

Безголовая версия здесь.

Когда дан вход 1:

вход-1

Когда дан вход 2:

вход-2

MD XF
источник
1
Это еще раз показывает превосходство Бейсика ....
Можно сохранить еще несколько байтов, установив для одной или нескольких переменных какое-либо часто используемое значение, например 10. Кроме того, предложите заменить его EXP(X)/(2*3.14159*S1*S1)наEXP(X)*M
потолок