برچسب: آمار

من همیشه فکر می کردم بقال محل از کجا می فهمه من مهندس هستم که بهم می گه «آقای مهندس چه شیری بدم خدمتون» و از اونطرف هم برام عجیب بود که یک آگهی «استخدام مهندس» با چند ده رزومه مواجه می شه که عملا اکثرشون توانایی‌های لازم رو ندارن. امروز فوربس در مطلب کشورهایی با بیشترین فارغ التحصیل مهندسی که عملا فقط یک اینفوگرافی بر اساس داده‌های فروم اقتصاد جهانی است جوابم رو داده: ایران از نظر تعداد فارغ التحصیل مهندسی در رتبه سوم جهان است! بعد از روسیه و آمریکا و البته حضور کشورهایی مثل اندونزی و ویتنام هم در بین ده کشور اول تعداد مهندس‌هایی که سالانه مدرکشون رو می گیرن، قابل توجه است.

حکایتی هست در این باره که راکفلکر شروع کارش از یک کارخونه کبریت سازی بوده. این داستان مدعی است که این مولتی میلیاردر که زمان جوانی در کارخونه کبریت سازی کار می کرده، روزی پیش مدیرش می ره و می گه ایده ای داره برای چند برابر کردن درآمد کارخونه و در مقابل این قول که در بخشی از این درآمد سهیم بده ایده اش رو می ده:

جعبه کبریت یک کشویی داره که کبریت ها اون تو هستن. این رو برعکس توی جعبه بذارین. هر کس جعبه رو باز می کنه بیشتر کبریت ها به زمین می ریزن و طرف خیلی زودتر از حالتی که کشویی درست نصب شده باشه، مجبور می شه یک جعبه کبریت دیگه بخره. چون کبریت چیز ارزونی است کسی به این امر توجه نخواهد کرد و هیچ کس هم شاکی نخواهد شد و شما چند برابر بیشتر جعبه کبریت خواهید فروخت.

ایده خبیثانه جالبیه ولی چرا اصولا از اول کبریت های کمتری توی جعبه نذاریم؟ اصولا اگر می شه اینجوری پول درآورد مطمئنا کشور ما جای بسیار مناسبی براشه.

برای بررسی این مساله به مغازه می ریم و یک بسته ده تایی کبریت توکلی که ظاهرا تنها بازیگر بزرگ کبریت در ایرانه رو می خریم. اونو به خونه می‌یاریم. ملحفه سبزمون رو رو زمین پهن می کنیم و به عنوان یک تنبل واقعی، فقط قسمتی که قراره تو عکس باشه رو اتو می کنیم:

و بعد می ریم سراغ ده تا جعبه کبریتی که خریدیم و پشتشون مدعی است که تعداد متوسط کبریت‌های هر جعبه چهل تا است:

هر جعبه کبریت رو جدا جدا روی قسمت اتو شده و صاف خالی کرده، عکس می گیریم:

اینکار رو برای هر ده تا جعبه کبریت انجام می دیم و اتاق پر از کبریت رو به قصد نشستن پشت کامپیوتر ترک میکنیم.

اول نیازمند یک برنامه هستیم که بتونه در یک عکس کبریت ها رو بشمره. بعد با یک اسکریپت این برنامه رو برای همه عکس ها اجرا خواهیم کرد. من برای اینکار اول کتابخونه‌های cv2 رو امتحان کردم ولی با numpy نتیجه بهتر و قشنگ تری گرفتم. در قدم اول کافیه عکس رو بخونیم، اون ۱ به کتابخونه می گه عکس رو بعد از خوندن خاکستری کنه که کار ما رو راحتتر می کنه:

kebrit = scipy.misc.imread(fileName, 1) # gray-scale image

و خروجی چیز شبیه این است:

اینکه تصویر رنگی است به خاطر اینه که به شکل پیش فرض از پالت jet برای نمایش استفاده کرده ایم که تصویر رو به صورت حرارتی نشون می ده. می گن این برای کارهای علمی راحتتر از خاکستری است. به راحتی می بینین که چیزی که می خونیم با اینکه در سبز معمولی یا خاکستری به نظر یکدست می‌یاد ولی پر از تیزی و نقطه نقطه است. مثل هر عکس دوربینی دیگه. برای حل این مشکل تصویر رو اسموت می کنیم. اسموت کردن در پردازش تصویر معمولا با استفاده از فیلتر گاوس به دست می یاد که کانولوشن گرفتنی است بین این ورودی و تابع گاوس. البته نگران نباشین، من چون سه بار معادلات دیفرانسیل رو افتادم و آخرش هم یکی دیگه به جام امتحان داد، پایه ام قوی شده (:

kebrit_smooth = ndimage.gaussian_filter(kebrit, 6)

حالا متغیر kebrit_smooth یک تصویر نرم شده از تصویر اصلی است:

قدم بعدی اینه که من بیام و یک ترش‌هولد روی ورودی بذارم و بگم بیخیال هر چیزی بشه که از این عدد خاص کوچکتر است. این عدد با کمی تجربه و کمی سعی خطا به دست اومده. دستور طبیعی باید این می بود:

tresh = 120
labeled, objectsNum = ndimage.label(kebrit_smooth < tresh)

تابع ndimage.label ورودی خودش رو بررسی می کنه و تعداد اجسام به هم پیوسته توی اون رو می شمره. اما مشکل اینه که حتی تک پیکسل های منفرد هم شمرده می شن. برای جلوگیری از این امر من بهش یک ساختار می دم که شکلم حداقل لازمه چنین ساختاری داشته باشه:

tresh = 120
removeOnes = np.ones((3,3), dtype="bool8")
labeled, objectsNum = ndimage.label(kebrit_smooth < tresh, structure=removeOnes)

دقت کنین که اون استراکچر در اصل می گه حداقل تصویر قابل تشخیص من باید این باشه:

array([[ True,  True,  True],
       [ True,  True,  True],
       [ True,  True,  True]], dtype=bool)

و نتیجه اش این می شه که نقطه های منفرد دیگه شمرده نشن. البته انتظارم اینه که با اسموت کردن مرحله قبل این نقطه ها به حداقل رسیده باشن. نتیجه این ترکیب چنین تصویری است:

و خب چنین خروجی متنی ای:

10boxes/IMG_2608.JPG , 29

حالا کافیه یک اسکریپت ساده، همه ده تا فایل رو به اون برنامه بده:

#!/bin/bash

for f in *JPG
do
    ../count.py $f
done

و به خروجی زیر برسیم:

$ ./doall.sh 
IMG_2606.JPG , 35
IMG_2607.JPG , 30
IMG_2608.JPG , 29
IMG_2609.JPG , 32
IMG_2610.JPG , 29
IMG_2611.JPG , 33
IMG_2612.JPG , 36
IMG_2613.JPG , 38
IMG_2614.JPG , 37
IMG_2615.JPG , 38

وقتشه سراغ زبون مورد علاقمون R بریم. یک زبان تخصصی برای کارهای آماری و وررفتن با اعداد و ماتریس ها و رفیق رفقاشون. فایل رو می خونیم و بخش مورد نظر رو جدا می کنیم:

> tavakoli <- read.csv(file="results.csv",head=FALSE,sep=",")
> matches <- tavakoli[2]$V2

حالا اطلاعات هر جعبه کبریت رو دارم. بذارین یک نگاه سریع بهش بندازیم:

> summary (matches)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  29.00   30.50   34.00   33.70   36.75   38.00 

بعله. می بینیم که میزان متوسط پایینتر از ۳۴ است. یعنی تعداد متوسط این ده قوطی کبریت ۶ تا کمتر از چیزی است که ادعا شده یا به بیان صحیح‌تر ۱۵٪ کمتر به ما جنس فروختن. جالبیش اینه حتی یک جعبه از این ده جعبه هم به چیزی که ادعا می شد میانگین است نرسیده. بذارین جلوتر بریم:

> hist(matches, xlim=c(28,42), ylim=c(0,4), main="هیستوگرام تعداد کبریت های توکلی در هر جعبه", sub="www.jadi.net", breaks=8)
> curve(dnorm(x, mean=mean(matches), sd=sd(matches))*10, add=TRUE, col="red", lwd=2) 

و روی این نمودار با کمی آمار احتمالات می تونیم احتمال اینکه در یک جعبه کبریت اتفاقی توکلی که خریده‌ایم و مدعی داشتن ۴۰ کبریت است، چهل یا بیشتر کبریت وجود داشته باشد را حساب کنیم:

> 1-pnorm(40, mean=mean(matches), sd=sd(matches)) #یک منهای سمت چپ نمودار نرمال در نقطه ۴۰ کبریت
[1] 0.03970968

بله. با اینکه کبریت توکلی مدعی است در هر جعبه اش تقریبا ۴۰ کبریت وجود داره، برنامه پردازش تصویر و نمودارهای نرمال ما نشون می دن که احتمال اینکه واقعا در یک جعبه کبریت توکلی چهل یا بیشتر کبریت باشه، سه صدم درصد بیشتر نیست

نتیجه‌ها

۰- موسسه استاندارد و حقوق مصرف کننده خاصی نداریم یا هنوز این مساله رو ندیدن
۱- راکفلر باید بیاد پیش تولید کننده‌های ما لنگ بندازه
۲- آمار شیرین و فان است
۳- می تونیم با کارهای علمی بامزه هم تفریح کنیم هم چیز یاد بگیریم

مرتبط

• هفته وایبر: منحنی نرمال، اسیدپاشی و اینکه چرا گروه‌های وایبری ترسناک می شن
• پست مهمان: رابطه استفاده از کلمه مهمل توی وبلاگ جادی و میزان بارون در تهران

و البته مثل همیشه هر جاییش ممکنه اشتباه داشته باشه. خوشحال می شم دوستان حرفه ای اصلاح کنن یا توسعه بدن و به این فکر کنیم که کاش آمار رو اینطوری به ما درس می دادن (:

برای من اینطرف کشیدن ساعت معنی اش اینه که دیگه نمی تونم لای پرده رو باز بذارم و به موقع با نور آفتاب بیدار بشم. الان نور افتاب منو ساعت شش بیدار خواهد کرد. اصولا این ساعت رو اینطرف و اونطرف کشیدن چه تاثیری توی نور روز داره و در سطحی بالاتر، من که این روزها از هفت از خونه بیرون می رم در بهترین حالت شش و ربع می رسم، چقدر از سال تو تاریکی می یام بیرون و تو تاریکی بر می گردم. بذارین به پیروی از این مقاله که کار مشابهی کرده با بسته نرم افزاری جادویی R، با اعداد ور بریم.

مشخصه که نیاز به کتابخونه ای داریم که فانکشن های sunrise و غیره رو دارن پس:

> library(maptools)

و کافیه لت و لانگ تهران رو بهش بدیم:

> portsmouth <- matrix(c(51.4231, 35.6961), nrow=1)

و تاریخ:

> for_date <- as.POSIXct("2014-09-25", tz="Asia/Tehran")

و ازش بخوایم ساعت طلوع خورشید امروز رو بده:

> sunriset(portsmouth, for_date, direction="sunrise", POSIXct.out=TRUE)
        day_frac                time
newlon 0.2460366 2014-09-25 05:54:17

بله ظاهرا درسته ! خورشید ما رو ساعت ۵ و ۵۴ دقیقه صبح بیدار می کنه که مثل گرگ پاشیم بریم سر کار. این کتابخونه همچنین می تونه به ما زمان غروب، طول روز و ظهر رو هم بده و از اون باحالتر می تونه ورودی هایی به فرم بردارهایی از روزها قبول کنه. پس این تابع:

# adapted from http://r.789695.n4.nabble.com/maptools-sunrise-sunset-function-td874148.html
ephemeris <- function(lat, lon, date, span=1, tz="UTC") {
 
  # convert to the format we need
  lon.lat <- matrix(c(lon, lat), nrow=1)
 
  # make our sequence - using noon gets us around daylight saving time issues
  day <- as.POSIXct(date, tz=tz)
  sequence <- seq(from=day, length.out=span , by="days")
 
  # get our data
  sunrise <- sunriset(lon.lat, sequence, direction="sunrise", POSIXct.out=TRUE)
  sunset <- sunriset(lon.lat, sequence, direction="sunset", POSIXct.out=TRUE)
  solar_noon <- solarnoon(lon.lat, sequence, POSIXct.out=TRUE)
 
  # build a data frame from the vectors
  data.frame(date=as.Date(sunrise$time),
             sunrise=as.numeric(format(sunrise$time, "%H%M")),
             solarnoon=as.numeric(format(solar_noon$time, "%H%M")),
             sunset=as.numeric(format(sunset$time, "%H%M")),
             day_length=as.numeric(sunset$time-sunrise$time))
 
}

می تونه لت و لانگ (طول و عرض) یک مکان رو بگیره و به اندازه date روز بعد از تاریخی که بهش گفتیم، به ما بگه که وضعیت نور روز چطوریه (اونم توی تایم زون مورد نظر):

> ephemeris(35.6961, 51.4231, "2014-09-25", 10, tz="Asia/Tehran")
         date sunrise solarnoon sunset day_length
1  2014-09-25     554      1156   1757   12.05001
2  2014-09-26     555      1155   1755   12.01276
3  2014-09-27     555      1155   1754   11.97550
4  2014-09-28     556      1155   1752   11.93826
5  2014-09-29     557      1154   1751   11.90104
6  2014-09-30     558      1154   1750   11.86384
7  2014-10-01     558      1154   1748   11.82666
8  2014-10-02     559      1153   1747   11.78952
9  2014-10-03     600      1153   1745   11.75241
10 2014-10-04     601      1153   1744   11.71535

جالبه ولی نه به اندازه کافی. مغز اکثر ما نمودارها رو بهتر از جدول ها می فهمه. پس ما با توابع ggplot دو تا نمودار می کشیم. یکی نمودار نواری که نشون بده چه بخشی از روز نور آفتاب داره و چه بخش هایی تاریکه و یک نمودار دیگه که به سادگی بگه هر روز از سال چند ساعت نور داره.

library(ggplot2)
library(scales)
library(gridExtra)

# create two formatter functions for the x-axis display
 
# for graph #1 y-axis
time_format <- function(hrmn) substr(sprintf("%04d", hrmn),1,2)
 
# for graph #2 y-axis
pad5 <- function(num) sprintf("%2d", num)
 
daylight <- function(lat, lon, place, start_date, span=2, tz="UTC", 
                     show_solar_noon=TRUE, show_now=TRUE, plot=TRUE) {
 
  stopifnot(span>=2) # really doesn't make much sense to plot 1 value
 
  srss <- ephemeris(lat, lon, start_date, span, tz)
 
  x_label = ""
 
  gg <- ggplot(srss, aes(x=date))
  gg <- gg + geom_ribbon(aes(ymin=sunrise, ymax=sunset), fill="#ffeda0")
 
  if (show_solar_noon) gg <- gg + geom_line(aes(y=solarnoon), color="#fd8d3c")
 
  if (show_now) {
    gg <- gg + geom_vline(xintercept=as.numeric(as.Date(Sys.time())), color="#800026", linetype="longdash", size=0.25)
    x_label = sprintf("زمان در لحظه ترسیم نمودار: %s", format(Sys.time(), "%Y-%m-%d / %H:%M"))
  }
 
  gg <- gg + geom_hline(yintercept=as.numeric("0700"), color="#1000F6", linetype="longdash", size=0.25)
  gg <- gg + geom_hline(yintercept=as.numeric("1815"), color="#1000F6", linetype="longdash", size=0.25)

  gg <- gg + scale_x_date(expand=c(0,0), labels=date_format("%b "))
  gg <- gg + scale_y_continuous(labels=time_format, limits=c(0,2400), breaks=seq(0, 2400, 200), expand=c(0,0))
  gg <- gg + labs(x=x_label, y="",
                  title=sprintf("طلوع و غروب خورشید در %s\n%s ", place, paste0(range(srss$date), sep=" ", collapse="تا ")))
  gg <- gg + theme_bw()
  gg <- gg + theme(panel.background=element_rect(fill="#525252"))
  gg <- gg + theme(panel.grid=element_blank())
 
  gg1 <- ggplot(srss, aes(x=date, y=day_length))
  gg1 <- gg1 + geom_area(fill="#ffeda0")
  gg1 <- gg1 + geom_line(color="#525252")
 
  if (show_now) gg1 <- gg1 + geom_vline(xintercept=as.numeric(as.Date(Sys.time())), color="#800026", linetype="longdash", size=0.25)
 
  gg1 <- gg1 + scale_x_date(expand=c(0,0), labels=date_format("%b "))
  gg1 <- gg1 + scale_y_continuous(labels=pad5, limits=c(0,24), expand=c(0,0))
  gg1 <- gg1 + labs(x="", y="", title="طول روز به ساعت")
  gg1 <- gg1 + theme_bw()
 
  if (plot) grid.arrange(gg, gg1, nrow=2)
 
  arrangeGrob(gg, gg1, nrow=2)
 
}

و حالا می تونیم هر چیزی بخوایم بکشیم.. مثلا برای تهران:

daylight(35.6961, 51.4231, "Tehran", "2014-03-21", 365, tz="Asia/Tehran")

باعث می شه به این برسیم:

نمودار بالا نشون می ده چه ساعت های از روز صبح است و چه ساعت هایی شب و نمودار پایینی نشون دهنده ساعت های روشنایی است که در هر روز داریم و توش به راحتی می تونین تاثیر اینطرف اونطرف کشیدن ساعت رو ببینین.

خط عمودی در نمودار بالا مشخص کرده الان کجای سال هستیم و دو خط نقطه چین افقی می گن اگر یک نفر لازم باشه از هفت صبح از خونه بره بیرون و شش و ربع برگرده، چه زمان‌هایی از نور آفتاب در زندگی اجتماعی دور است. این نمودار ممکنه باعث بشه شما از کارتون استعفا بدین.

می بنین اعداد چقدر قشنگن؟

نقطه اشتراک سوریه، ترکمنستان، کره شمالی و اریتره اینه که در گزارش امسال وضعیت آزادی بیان در جهان، از بین ۱۷۹ کشور کره زمین، تنها چهار کشوری هستن که از ایران پایین تر هستن. به عبارت دیگه چیزی که این چهار اسم رو به هم پیوند می ده اینه که تنها کشورهای جهان هستن که توشون آزادی بیان بیشتر از ایران تحت فشار است.